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1  Summary 

This  research  considers  the  mathematical  representation  of  complex  network  environments 
to  enable  network  reconstruction  and  analysis.  In  this  study,  a  complex  network  environment 
is  the  interconnection  of  controlled  dynamical  systems  resulting  in  both  nonlinear  and  noisy 
behavior  as  well  as  complicated  interconnection  structure.  Classical  representations  of  such 
systems,  such  as  coupled  differential  equations,  are  not  effective  for  our  purposes  because  they 
are  too  detailed,  resulting  in  a  heavy  information  load  on  the  experimental  data  required  for 
reconstruction.  Moreover,  the  typical  simplification  found  by  grouping  parts  of  the  system 
into  well- characterized  subsystems  and  then  restricting  the  desired  structural  information  to 
be  only  the  interconnection  structure  among  subsystems  is  equally  problematic,  since  even 
here  one  must  determine  how  to  partition  all  system  states  (even  the  unmeasured  ones)  into 
the  right  subsystems;  this  can  be  a  tremendous  informational  burden  or  even  impossible. 

Here,  we  leverage  a  new  idea  about  how  to  represent  system  structure  to  reduce  the 
information  burden  on  the  experimental  data  required  for  network  reconstruction.  In  this 
new  representation,  we  look  for  the  open-loop  causal  dependencies  among  measured  vari¬ 
ables.  As  a  result,  distinct  dependencies  between  different  variables  may  actually  depend  on 
the  same  underlying  unmeasured  variable  and  thus  be  entangled.  In  this  situation,  this  new 
representation  would  not  distinguish  between  entangled  components  and  non-entangled  com¬ 
ponents.  Nevertheless,  the  benefit  of  remaining  agnostic  about  the  potential  entanglement 
of  system  dynamics  through  unmeasured  variables  is  that  this  new  type  of  structure  can  be 
derived  from  O(n)  experiments,  where  n  is  the  number  of  measured  variables.  Algorithms 
for  accomplishing  this  reconstruction  are  discussed,  with  special  emphasis  on  how  to  han¬ 
dle  noise  and  underlying  nonlinearities.  Also,  further  theoretical  analysis  reveals  interesting 
characteristics  about  minimal  representations  of  these  new  models  and  their  potential  for 
conducting  structure-dependent  robustness  analyses.  This  robustness  analysis  provides  the 
basis  for  exploring  vulnerability  and  network  security. 

A  unique  strength  of  this  work  is  its  equal  effort  on  both  theoretical  development  and  ex¬ 
perimental  validation.  Using  wireless  mesh  networks  as  an  application  testbed,  the  “network 
design  cycle”  is  employed  to  iteratively  develop  theoretical  models,  rate  control  protocols, 
and  experimental  validation  in  a  process  that  lays  the  foundation  for  the  development  of 
a  model  of  wireless  mesh  networks  suitable  for  the  new  network  reconstruction  methods. 
Although  not  yet  complete,  this  work  launches  new  possibilities  for  wireless  mesh  networks 
where  system  intruders  can  be  automatically  detected  by  the  way  they  induce  coupling 
among  link  rates  in  the  network.  This  “packet-free”  approach  is  distinct  from  current  tech¬ 
niques  and  suggests  a  bold  new  research  program.  Bio-chemical  reaction  networks  provide 
another  testbed  for  theoretical  validation,  and  that  work  is  just  emerging. 

2  Introduction 

The  study  of  network  environments  have  received  considerable  attention  from  various  disci¬ 
plines  over  the  last  20  years.  This  is  in  part  due  to  the  rise  of  the  internet,  but  it  also  is  due 
to  advances  in  biology  that  have  brought  new  kinds  of  networked  “machines”  into  clear  view 
that  contrast  sharply  with  our  solid-state,  engineered  systems.  The  research  discussed  here 
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takes  a  fresh  view  of  these  new,  “complex”  network  environments  and  answers  fundamental 
questions  about  1)  how  to  model  them,  2)  the  design  of  experiments  necessary  to  discover 
their  structure  (and  thus  adapt  system  inputs  to  optimize  the  resulting  performance),  and 
3)  the  relationship  between  network  structure  to  vulnerability  and  attack. 

Specifically,  this  work  explores  these  issues  in  the  context  of  both  wireless  mesh  and  bio¬ 
chemical  reaction  networks.  Although  wildly  different  application  areas,  this  research  unites 
theoretical  work  that  clarifies  fundamental  limitations  of  complex  networks  with  network 
engineering  and  systems  biology  to  implement  specific  designs  and  experimentally  verify  the 
theoretical  discoveries  in  what  we  call  the  “network  design  cycle.” 

2.1  What  is  a  “complex”  networked  environment? 

The  descriptor  “complex”  has  been  used  in  other  studies  to  characterize  networks  that  are 
large,  meaning  that  a  graph  used  to  represent  the  network  has  a  large  number  of  nodes. 
Often  the  number  of  edges  is  also  taken  to  be  large  and  devoid  of  certain  regular  structure, 
so  that,  for  example,  a  tree-structured  graph  would  not  be  considered  “complex.”  In  these 
studies  one  typically  looks  for  hidden  regularities  or  patterns  that  characterize  the  structure 
of  the  graph  in  simple  terms  in  spite  of  these  other  “complexities.”  A  typical  example  would 
be  small-world  networks. 

In  this  work,  we  consider  different  environments,  where  “network”  refers  to  an  inter¬ 
connected  dynamical  system.  In  these  environments,  system  behavior  is  as  important  a 
characteristic  of  the  network  as  is  its  structure.  The  descriptor  “complex”  then  refers  to 
both  the  network  behavior  as  well  as  its  structure. 

This  work  addresses  networks  with  complex  behavior  by  considering  systems  with  un¬ 
derlying  nonlinear  dynamics  and  noisy  measurements.  As  a  first  step,  the  results  presented 
here  leverage  Lyapunov  theory  to  apply  linear  analysis  locally  to  regions  near  equilibria  of 
the  underlying  nonlinear  system.  Extending  these  results  globally  to  the  underlying  non¬ 
linear  system  is  not  immediately  obvious  and  requires  new  thinking  about  the  meaning  of 
structure,  beyond  that  already  developed  in  this  research.  As  a  result,  global  analysis  is  left 
for  future  research. 

This  work  also  addresses  networks  with  complex  structure.  Besides  the  usual  definition, 
of  a  “complex”  network  structure  represented  by  a  graph  with  a  large  number  of  nodes  and 
non-trivial  edge  patterns  (e.g.  allowing  for  arbitrarily  complicated  feedback  relationships), 
this  work  also  makes  a  particular  contribution  by  developing  representation  and  analysis  tools 
applicable  to  networks  of  dynamical  systems  with  potential  hidden  entanglements  among  un¬ 
measured  variables.  This  “potential  entanglement”  type  of  network  complexity  is  previously 
unaddressed  in  the  literature,  yet  it  becomes  particularly  important  for  inferring  network 
structure  from  behavioral  data. 

Appreciating  the  power  of  structural  representations  that  allow  for  potential  entangle¬ 
ment  among  unmeasured  variables  to  simplify  network  inference  problems  is  subtle,  but  it 
is  a  central  contribution  of  this  research  program.  Consider,  for  example,  a  network  that 
is  composed  of  the  interconnection  of  various  subsystems.  By  definition,  each  subsystem’s 
internal  states  only  affect  that  subsystem,  and  the  interconnection  variables  are  themselves 
measured  quantities.  Inferring  this  network  “subsystem”  interconnection  structure  from  data 
thus  demands  the  discovery  of  the  true  partition  of  all  of  the  system’s  unmeasured  states 
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into  their  appropriate  subsystem  components.  Discovering  such  a  partition  from  measured 
data  can  be  extremely  difficult  or  impossible.  This  research,  on  the  other  hand,  leverages 
a  different  representation  of  system  structure,  called  signal  structure,  that  does  not  rely  on 
the  idea  of  subsystems  and  allows  for  potential  entanglement  among  unmeasured  states.  As 
a  result,  inferring  a  system’s  signal  structure  requires  much  less  information,  and  thus  fewer 
experiments,  than  inferring  a  system’s  subsystem  structure. 

One  contribution  of  this  research  is  to  thoroughly  understand  the  relationships  between 
a  system’s  subsystem  and  signal  structures.  Often,  systems  with  solid-state  components 
(such  as  routers  in  the  internet)  have  subsystem  and  signal  structures  that  are  equivalent. 
Sometimes,  however,  systems  have  a  fluid-like  character  (such  as  bio-chemical  reaction  net¬ 
works  or  wireless  mesh  networks)  and  the  resulting  subsystem  and  signal  structures  can  be 
very  different.  Characterizing  informativity  conditions  and  developing  scalable  algorithms 
for  network  reconstruction  (i.e.  inference)  of  signal  structure  are  the  primary  theoretical 
contributions  of  this  effort.  These  models  of  complex  networked  environments  also  facilitate 
a  novel  robustness  analysis  that  leads  to  new  results  about  system  vulnerability  and  secu¬ 
rity.  These  contributions  are  highlighted  when  applied  to  complex  network  environments 
that  exhibit  both  the  behavioral  and  structural  complexity  as  described. 

2.2  What  is  the  network  design  cycle? 

Besides  the  theoretical  contributions  of  this  work,  this  research  represents  an  active  col¬ 
laboration  between  theoretical  development  and  physical  implementation  and  testing  on 
real  complex  networks.  This  collaboration  is  most  evident  in  our  work  on  wireless  mesh 
networks,  where  active  modeling  and  protocol  design  efforts  lead  to  simulation,  implementa¬ 
tion,  and  experimental  testing  on  our  live  wireless  mesh  testbed.  Similar  collaborations  for 
bio-chemical  reaction  networks  began  to  emerge  during  this  study,  but  the  implementation 
and  experimental  testing  is  incomplete  and  part  of  ongoing  research. 

The  network  design  cycle  defines  the  scientific  process  we  engage  that  unites  our  theoret¬ 
ical  and  applied  work.  The  next  section  discusses  each  part  of  the  design  cycle  in  detail  as 
part  of  our  research  methods,  assumptions,  and  procedures.  Section  7  the  offers  highlights 
of  the  research  program  results,  beginning  with  specific  applications  and  ending  with  general 
theoretical  results.  Sections  7.1  and  7.2  detail  our  development  of  new  models  for  wireless 
mesh  networks,  leading  to  new  rate  control  protocols  that  can  drastically  improve  network 
performance.  Section  7.3  and  7.4  then  demonstrate  our  network  reconstruction  efforts  ap¬ 
plied  to  nonlinear  bio-chemical  reaction  networks  with  noisy  measurements.  Section  7.5  then 
highlights  the  general  theoretical  underpinnings  regarding  the  meaning  of  structure  in  inter¬ 
connected  dynamical  systems  and  the  mathematical  relationships  among  different  types  of 
structure.  Section  8  then  concludes  the  report. 

3  Methods,  Assumptions,  and  Procedures 

To  address  these  problems,  our  work  uses  a  method  we  call  the  Network  Design  Cycle,  as 
shown  in  Figure  1.  We  start  by  formulating  a  mathematical  model  of  the  network,  precisely 
characterizing  how  it  operates.  In  a  wireless  network,  this  may  involve  describing  how  nodes 
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Figure  1:  The  Network  Design  Cycle 


interact  via  carrier  sensing  and  interference.  Based  on  this  model,  we  carefully  formulate 
the  research  problem,  such  as  providing  optimal  performance  or  detecting  malicious  nodes. 
The  next  step  is  to  design  algorithms  and  protocols  that  will  solve  the  problem.  We  then 
implement  the  solution,  which  often  requires  considering  additional  complexities  that  arise 
when  building  and  deploying  code  on  a  network.  Finally,  we  run  experiments  to  validate  that 
the  solution  works  as  designed,  for  example  by  providing  optimal  performance  or  security 
guarantees.  If  performance  deviates  from  what  is  expected,  this  means  that  our  original 
model  or  problem  formulation  was  wrong.  This  leads  to  another  iteration  of  the  design 
cycle,  where  we  create  a  refined  model  or  reformulate  the  problem. 

For  the  steps  that  involve  implementation  and  experimental  evaluation,  we  use  a  variety 
of  tools.  We  use  MATLAB  to  numerically  evaluate  algorithms,  ensuring  that  they  converge 
to  the  expected  solution  and,  in  some  cases,  provide  optimal  performance.  For  experimental 
testing,  we  use  a  wireless  mesh  network  deployed  in  our  department’s  building,  consisting 
of  30+  computers  configured  as  wireless  routers,  spread  over  two  floors.  The  mesh  network 
can  be  configured  so  that  it  uses  802.11a,  and  we  ensure  that  we  select  a  channel  with  no 
other  traffic,  so  that  we  can  run  experiments  without  any  outside  noise  when  desired.  Our 
protocols  are  implemented  using  WiFu  [4],  a  toolkit  developed  with  funding  from  NSF  that 
enables  rapid  deployment  of  experimental  wireless  protocols  in  user-space. 

4  Results  and  Discussion 

4.1  Modeling  Wireless  Networks  with  Partial  Interference 

Because  wireless  networks  use  shared  communication  channels,  contention  and  interference 
can  significantly  degrade  throughput  and  fairness.  A  common  approach  to  solving  this 
problem  is  to  form  a  model  representing  the  wireless  network  that  imposes  constraints  on 
the  rates,  and  then  design  the  controller  to  compute  a  distributed  algorithm  that  solves 
an  optimization  problem.  In  wireless  networks,  these  constraints  are  primarily  imposed  by 
the  channel  capacity  among  competing  nodes  [6].  By  formulating  the  problem  so  that  it  is 
convex,  well  established  techniques  can,  in  most  cases,  be  applied  to  design  a  distributed 
algorithm  that  computes  the  rates  [13]. 

Traditional  graph-based  models  impose  resource  constraints  using  a  contention  graph, 
where  communication  links  between  nodes  are  represented  by  vertices  and  contention  is 
represented  by  edges.  Two  links  are  said  to  contend  with  each  other  if  they  caimot  be  active 
at  the  same  time  without  causing  collisions.  Existing  graph-based  models  conservatively 
characterize  interference  as  a  binary  effect  [5,  10].  Under  this  model,  an  interfering  node  is 
assumed  to  corrupt  all  of  the  packets  received  at  a  remote  node,  while  non-interfering  nodes 
have  no  effect.  Binary  interference  is  represented  in  a  contention  graph  by  simply  treating 
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Figure  2:  A  network  topology  graph. 


Figure  3:  A  contention  graph  for  the  sample  network  topology,  based  on  the  partial  inter¬ 
ference  model. 

it  as  contention,  that  is,  if  one  link  interferes  with  another,  neither  may  send  at  the  same 
time. 

4.1.1  Model 

While  this  approach  does  lead  to  convex  problem  formulations,  it  does  not  accurately  rep¬ 
resent  the  interference  relationship  among  links.  Contention  occurs  when  sending  nodes  can 
sense  one  another  and  thus  take  turns,  whereas  interference  occurs  when  two  non-contending 
senders  transmit  simultaneously,  causing  a  packet  to  be  corrupted  at  a  receiving  node.  A  re¬ 
cent  measurement  study  reveals  that  interference  is  typical  and  partial  [14].  This  means  that 
transmissions  from  an  interfering  node  may  corrupt  only  a  fraction  of  the  packets  received 
at  a  remote  node.  Modeling  interference  as  contention  results  in  a  misleading  optimization 
problem,  where  capacity  may  be  wasted  and  actual  receiving  rates  may  be  far  from  fair. 

We  have  developed  a  new  model  of  the  wireless  mesh  network  that  accurately  models 
partial  interference.  In  this  model,  contention  is  still  represented  as  an  undirected  edge 
between  two  vertices  (links),  however  interference  is  modeled  as  a  directional  edge  from  the 
interfering  link  to  the  receiving  link  that  is  affected  by  the  interference.  Figure  2  shows 
a  simple  topology  of  a  wireless  mesh  network,  with  the  resulting  contention  graph  shown 
in  Figure  3.  From  this  graph,  we  find  the  set  of  maximal  cliques,  which  results  in  clique 
constraints.  In  our  example  of  Figure  3,  the  constraints  for  Clique  1  and  Chque  2,  are  the 
same  as  in  the  binary  interference  model,  but  the  constraint  between  links  2  and  3  is  modeled 
separately. 

Interference  relationships  impose  constraints  on  the  receiving  rates,  as  illustrated  in  the 
figure,  where  rt  is  the  effective  receiving  rate  of  link  /,  st  is  the  sending  rate  of  link  /,  dt  is  the 
inherent  loss  of  the  link  (e.g.  due  to  obstacles  or  noise),  Cj  is  the  capacity  of  chque  j,  and 
the  term  (1  —  o^s,)  is  the  loss  of  link  l  due  to  interference  from  an  interfering  node  i.  This 
constraint  is  taken  from  a  recent  measurement  study  of  wireless  mesh  networks  showing  that 
partial  interference  from  one  link  can  be  modeled  as  a  linear  function  [14].  The  interference 
factor  an  represents  the  degree  of  partial  interference  inflicted  by  the  interferer.  It  is  in  the 
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range  0  <  an  <  1,  and  is  directional,  meaning  that  an  may  be  significantly  different  from  an. 
Interfering  factors  can  be  experimentally  measured  between  any  pair  of  links  in  a  network 
by  methods  suggested  in  [14,  16],  constructing  an  interference  map  for  the  network.  A  study 
shows  that  interfering  transmissions  are  independent  of  each  other,  and  the  joint  impact  of 
interferers  to  a  receiving  node  is  merely  the  product  of  their  isolated  impacts  [14], 


4.1.2  Problem  Formulation 

Given  a  contention  graph  with  maximal  cliques  C  and  an  interference  map  A,  we  formulate 
an  optimization  problem  that  maximizes  the  sum  of  link  utilities,  which  are  functions  of  link 
receiving  rates,  in  a  wireless  mesh  network: 


P  :  max/(r)  =  V  IJ  (ri) 
leL 


(i) 


subject  to: 

Si  >0,  V/  G  L , 

ri  =  dtsi  (1  -  auSi),  VI  G  L, 

iei(i) 

Y  si-  cv  vi  e  c- 


(2) 

(3) 

(4) 


We  assume  the  utility  function  U  of  a  link  is  continuously  differentiable,  strictly  concave, 
monotonically  increasing,  and  approaches  negative  infinity  as  the  argument  approaches  zero 
from  the  right. 

Assuming  a  logarithmic  utility  function,  and  ignoring  the  ratio  di  because  it  does  not 
affect  the  optimality  of  our  solution,  we  can  formulate  this  as  a  convex  problem  P7 


P7  :  max/(s) 


(5) 


where 


subject  to: 


f\s)  =  El  In  si  +  Y  ln  (1  ~  ansi)  )  • 

ieF(l) 


(6) 


leL 


si  >  0,  V/  G  L,  (7) 

Y  Sl-  G’  e  c ■  (8) 

Note  that  a  similar  formulation  that  optimizes  flow  receiving  rates,  rather  than  link 
receiving  rates,  is  not  convex.  Finding  a  convex  reformulation  is  a  part  of  our  ongoing  work. 


APPROVED  FOR  PUBLIC  RELEASE:  DISTRIBUTION  UNLIMITED 

6 


4.1.3  Algorithm 

Because  the  link  formulation  of  this  problem  is  convex,  it  is  straightforward  to  derive  a 
distributed  algorithm  to  solve  it,  based  on  the  methods  presented  in  [13].  We  use  the  gradient 
projection  method  to  iteratively  obtain  the  optimal  price,  A,  for  the  problem.  Using  a  step 
size  7  in  the  negative  direction  of  the  gradient  gives  the  algorithm 

A j(k  +  1)  =  max  I  0,  A j(k)  -  7 (cj  -  ^  si(X(k)))  J  ,  (9) 

\  ieL(j)  / 

where 

si( A)  =  arg  max g(sh  A).  (10) 

Si 

and  _ 

g(sh  A)  =  In  si  +  ^  In  (1  -  ausi)  -  st  ^  A  j,  (11) 

ieF(i)  jeC(i) 

where  F(l)  is  the  set  of  all  links  that  interfere  with  /  and  C(l)  is  the  set  of  all  links  in  a 
clique  with  l. 

This  means  that  each  link  in  a  maximal  clique  can  share  with  each  other  their  current 
rates,  sj,  which  is  only  local  information,  and  then  compute  the  next  price  A  j,  which  in  turn 
leads  each  link  to  calculate  new  prices.  The  convergence  of  the  algorithm  is  well  established 
in  the  literature,  even  when  it  is  asynchronous.  Once  A  converges  to  the  optimal  solution, 
the  optimal  solution,  s*,  is  given  by 

s*  =  s(A*). 


4.1.4  Numerical  Results 

We  use  MATLAB  to  demonstrate  the  situations  in  which  the  partial  interference  (PI)  model 
outperforms  the  binary  interference  (BI)  model,  and  by  how  much.  For  binary  interference, 
we  assume  that  the  model  can  either  model  interference  as  contention  or  ignore  it,  depending 
on  which  gives  the  largest  utility.  In  most  cases,  the  BI  model  ignores  interference  when  it 
is  low  and  models  it  as  contention  when  it  is  high.  To  compare  the  models,  we  consider 
the  ratio  R  of  performance  between  the  PI  and  BI  models,  using  three  topologies,  each  with 
clique  capacities  of  0.85. 

Figure  4a  shows  the  first  topology,  where  /  links  interfere  with  a  single  link  with  a 
common  interference  factor  a,  but  do  not  interfere  with  each  other.  Figure  5a  shows  the 
performance  ratio  for  this  topology,  with  the  dotted  curve  showing  when  R  begins  to  be 
greater  than  one.  Interestingly,  the  PI  model  and  the  BI  model  perform  exactly  the  same 
for  values  of  a  below  0.59.  This  is  because,  for  low  values  of  a,  the  cost  of  interference  is 
offset  by  the  gain  of  the  interferer  sending  at  full  capacity.  Thus,  both  the  PI  model  and  the 
BI  model  calculate  sending  rates  at  full  capacity  for  each  link.  For  larger  values  of  a  and  /, 
the  PI  model  outperforms  BI  more  than  1.7  times. 

Figure  4b  shows  the  second  topology,  where  a  single  link  has  interference  factor  a  on  N 
links  that  contend  in  a  single  clique.  Figure  5b  plots  the  performance  ratio  for  this  topology, 
with  the  dotted  curve  showing  where  R  begins  to  be  greater  than  one.  The  PI  model  starts 
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Clique 


(a)  I  links  interfering  with  a 
single  link. 


Figure  4:  Topologies 


(b)  N  contenders  in  a  clique 
with  one  interferer. 

for  numerical  results. 


used 


(a)  Ratio  R  of  performance  be¬ 
tween  the  PI  model  and  the  BI 
model  for  I  intererers  on  one  link. 


(b)  Ratio  R  of  performance  be¬ 
tween  the  PI  model  and  the  BI 
model,  for  N  contenders  with  one 
interferer. 


(c)  Ratio  R  of  performance  be¬ 
tween  the  PI  model  and  the  BI 
model,  for  I  interferer s  and  N  con¬ 
tenders,  with  a  =  0.4. 


Figure  5:  Numerical  results  for  PI  model. 


performing  better  than  the  BI  model  at  much  lower  values  of  a  when  N  is  large.  This  is  due 
to  the  fact  that  the  contending  links  already  have  small  rates  as  a  consequence  of  sharing 
the  medium.  Utilities  are  lowered  much  more  by  interference  when  sending  rates  are  small. 
Thus,  even  for  low  values  of  a,  the  PI  model  does  not  calculate  sending  rates  at  full  capacity. 
However,  for  higher  values  of  a  and  N,  the  PI  model  outperforms  the  BI  model  only  about 
1.1  times. 

To  further  demonstrate  the  worth  of  the  PI  model,  we  consider  a  topology  combining 
features  of  the  first  two,  where  I  links  have  a  fixed  interference  factor  a  =  0.4  on  N  links 
that  contend  in  a  single  clique.  Figure  5c  plots  R  for  this  topology.  Experimental  results 
show  that  it  is  typical  for  interference  factors  to  range  anywhere  between  zero  and  one  in 
a  real  network,  with  usually  at  least  one  interferer  on  a  link  having  a  factor  of  at  least 
a  =  0.8,  so  choosing  a  =  0.4  in  this  topology  is  a  somewhat  conservative  comparison  [14]. 
The  combined  effect  of  several  interferers  and  several  contenders  causes  the  PI  model  to 
perform  significantly  better  than  the  BI  model. 

4.1.5  Implementation 

To  validate  the  performance  benefits  of  the  PI  model,  we  implement  a  rate  controller  as 
an  application  for  the  Linux  operating  system,  using  the  WiFu  toolkit  [4].  The  primary 
challenges  for  the  implementation  are  (1)  constructing  a  network  interference  map,  giving 
the  contention  and  interference  variables  we  use  for  our  partial  interference  model,  (2)  based 
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on  the  interference  map,  forming  a  contention  graph  and  enumerating  maximal  cliques, 
(3)  implementing  the  rate  controller  algorithms,  and  (4)  scheduling  flows  according  to  the 
calculated  rates. 

Interference  Map  To  measure  the  network  interference  map,  we  use  unicast  transmis¬ 
sions  between  pairs  of  nodes,  similar  to  the  method  discussed  by  Padhye  et  ah  [15].  The 
measurement  is  performed  in  four  steps: 

1.  Sending  nodes  select  receivers.  Each  sending  node  selects  a  receiver  to  form  two  sep¬ 
arate  unicast  links,  say  A— >C  and  B— >D.  The  receivers  should  be  chosen  in  a  way 
that  transmissions  from  one  sending  node  do  not  interfere  with  the  packet  reception 
of  another  link.  This  is  to  separate  the  impact  of  interference  from  contention. 

2.  Sending  nodes  take  turns  to  transmit  data  to  their  respective  receivers  using  unicast 
at  full  link  capacity.  Receivers  calculate  the  transmission  rate.  In  the  example,  C 
calculates  Ra  and  D  calculates  R/, . 

3.  Sending  nodes  concurrently  transmit  data  to  their  respective  receivers  using  unicast  at 
full  link  capacity.  Receivers  calculate  the  transmission  rate.  In  the  example,  C  calcu¬ 
lates  Rba  and  D  calculates  R The  node  in  the  superscript  is  the  potential  contending 
node. 

4.  Receivers  calculate  the  contention  ratio.  For  node  C,  the  contention  ratio  is  defined  as 

JDb  *  ffQ' 

s2-,  and  D  calculates  the  ratio  as  Note,  the  two  contention  ratios  could  be  different 

lta  rxb 

for  asymmetric  contion.  Node  B  contends  with  node  A  if  ^  <  1,  and  similarly  node 

rta 

Ra 

A  contends  with  B  if  <  1. 

rib 

The  above  steps  are  repeated  for  each  permutation  of  nodes  of  interest  for  contention 
measurement.  A  pair  of  nodes  is  considered  to  be  contending  with  each  other  if  either  one  of 
the  contention  ratios  is  less  than  1,  because  the  partial  interference  model  assumes  symmetric 
contention  between  a  pair  of  neighboring  nodes. 

The  interference  map  is  measured  using  a  similar  unicast  approach.  However,  the  inter¬ 
fering  node  should  be  within  the  interference  range  of  the  receiver  and  outside  of  the  carrier 
sense  range  of  the  sender  of  the  interferee  link.  To  measure  interference  of  node  B  to  link 
A— )-C,  node  A  Erst  transmits  to  C  at  link  capacity  while  B  remains  silent,  and  C  calculates 
the  receiving  rate  from  A,  Ra.  In  the  next  step,  both  A— >-C  and  B— )-D  transmit  at  link 
capacity  concurrently,  and  C  calculates  the  receiving  rate  Rba  when  transmissions  from  B  are 
present.  Node  B  interferes  with  link  A— >-C  if  ^  <  1,  and  the  interference  factor  when  B 

transmits  at  full  link  capacity  is  defined  as  ^1  —  . 

Enumerating  Maximal  Cliques  To  calculate  fair  rates,  links  need  to  construct  their 
local  contention  graphs  and  hnd  the  maximal  cliques  using  the  network  interference  and 
contention  map.  Unfortunately,  enumerating  maximal  cliques  in  an  arbitrary  graph  is  a 
well-known  NP-hard  problem,  and  the  problem  is  extremely  difficult  in  dense  graphs.  With 
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the  binary  interference  model,  links  within  interference  range  of  each  other  cannot  be  active 
concurrently.  As  a  result,  the  contention  graph  is  likely  to  be  dense,  because  there  are 
likely  a  large  number  of  interferers  to  a  remote  node.  Existing  graph-based  proposals  adopt 
approximations  by  either  over-conservatively  assigning  links  within  2  hops  apart  to  the  same 
maximal  clique  [2],  or  requiring  major  modifications  to  underlying  link  layer  protocols  [9,  5]. 

With  the  partial  interference  model,  it  is  viable  to  design  efficient  and  accurate  protocols 
that  enumerate  maximal  the  cliques  in  a  contention  graph.  The  work  presented  in  [7]  suggests 
that  efficient  algorithms  exist  for  enumerating  maximal  cliques  in  graphs  that  are  1)  sparse 
and  2)  closed  under  the  operation  of  taking  subgraphs.  The  maximal  cliques  in  the  partial 
interference  model  only  consist  of  links  with  senders  within  the  carrier  sense  range  of  each 
other.  In  a  typical  wireless  mesh  network,  the  number  of  immediate  contenders  is  quite 
limited,  and  we  can  assume  the  contention  graph  satisfies  the  above  two  assumptions. 

Accordingly,  we  use  the  Bron-Kerbosch  algorithm  [3]  to  calculate  all  the  maximal  cliques 
that  a  link  belongs  to.  The  algorithm  is  executed  within  each  node,  and  the  CPU  overhead 
is  trivial  as  compared  to  the  wireless  communication  overhead.  Our  method  is  able  to  work 
in  real  time  as  compared  to  previous  work  in  the  literature  that  conservatively  assigns  any 
node  within  two  hops  to  the  same  clique. 

Rate  Controller  Algorithms  When  implementing  the  rate  controller  algorithms,  there 
are  a  number  practical  issues  we  must  consider.  First,  the  link-based  formulation  of  the 
partial  interference  model  doesn’t  capture  what  we  really  value  -  flow  utility  rather  than 
link  utility.  To  overcome  this  limitation,  we  set  a  weight  on  each  link  equal  to  the  number  of 
flows  that  traverse  the  link.  Links  with  more  flows  will  have  heavier  weights  in  the  optimal 
rate  control,  and  thus  be  assigned  more  bandwidth.  To  prevent  unfairness  between  flows 
sharing  the  same  link,  each  link  equally  divides  the  assigned  bandwidth  among  its  flows. 

Second,  the  partial  interference  model  is  based  on  the  assumption  that  each  link  has  an 
infinite  backlog  of  packets  to  send,  which  clearly  may  not  hold  in  practice.  This  means  that 
when  a  link  is  assigned  a  given  fair  rate,  it  may  not  actually  be  able  to  transmit  packets 
at  that  rate.  To  deal  with  this  problem,  our  implementation  uses  the  actual  transmission 
rates  for  each  link  when  calculating  prices,  rather  than  the  assigned  rate.  If  a  sending  node 
transmits  at  a  lower  rate  than  expected  due  to  insufficient  packets,  the  sum  rate  of  the  clique 
becomes  lower  than  the  optimal  target,  and  the  clique  price  becomes  lower,  leading  all  links 
in  the  clique  to  increase  their  rates.  In  this  way,  flows  with  sufficient  packets  are  able  to 
utilize  the  idle  channel.  However,  the  rate  allocation  is  unable  to  converge  to  the  optimal 
target,  because  there  is  always  a  gap  between  the  actual  and  the  optimal  transmission  rate 
because  of  the  insufficient  packets. 

Scheduling  Flows  Once  the  rate  controller  has  determined  the  optimal  rate  for  each  link, 
we  need  to  schedule  the  packets  for  each  flow  traversing  each  link,  ensuring  that  the  total 
rate  of  all  flows  on  a  link  obey  the  rate  limit.  To  do  this,  we  use  WiFu,  which  allows  user- 
space  programs  on  a  Linux  system  to  intercept,  modify,  and  reschedule  packets  as  they  are 
forwarded.  Figure  6  shows  the  architecture  for  our  system. 

WiFu  intercepts  packets  using  the  Linux  iptables  software  with  the  netfilter  interface 
[1],  The  intercepted  packets  are  then  given  to  a  handler,  depending  on  where  the  packet  was 
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WiFu  Toolkit 


Figure  6:  Fair  rate  control  implementation  architecture 


intercepted,  relative  to  the  forwarding  process.  Control  packets  are  send  to  a  separate  handler 
that  the  rate  controller  uses  to  calculate  fair  rates.  These  messages  include  notification  from 
neighbors  regarding  their  current  link  rates,  which  the  rate  controller  uses  to  calculate  a 
clique  price.  The  rate  controller  uses  a  per-flow  scheduler,  with  logical  queues  for  each 
flow  that  traverses  a  link.  Data  packets  are  placed  in  a  shared  physical  queue,  and  then 
transmitted  by  the  scheduler  according  to  their  assigned  rate. 

4.1.6  Experiments 

We  evaluate  the  partial  interference  model  on  a  wireless  mesh  testbed  located  in  our  de¬ 
partment’s  building.  One  part  of  this  testbed  is  shown  in  Figure  7.  We  run  a  number  of 
experiments,  but  in  this  report  focus  a  scenario  where  two  connections  are  active:  mesh25 
— >•  mesh  31  and  meshl7  — >  mesh27.  Due  to  the  topology  and  signal  strengths,  the  upper 
connection  in  this  figure  interferes  with  the  lower  connection,  with  an  interference  factor 
a  =  0.546.  This  means  that  at  full  power  the  upper  connection  corrupts  an  average  of  54.6% 
of  the  packets  being  sent  on  the  lower  connection. 

As  we  might  expect,  the  solution  to  the  rate  optimization  problem  in  this  case  requires 
having  the  interferer  slow  down  its  sending  rate  to  improve  the  performance  of  the  interferee 
and  maximize  the  sum  of  the  link  utilities.  In  this  case,  the  upper  connection  should  send 
at  a  normalized  rate  of  0.915  (1052  KBps)  while  the  lower  connection  maintains  a  rate  of 
1.0  (1150  KBps).  The  receiving  rate  of  the  upper  link  should  be  1052  Kbps  while  the  lower 
rate  should  receive  data  at  a  rate  of  575  Kbps.  Figure  8a  shows  the  actual  received  rates  for 
each  connection.  Note  that  the  lower  connection  receives  only  475  KBps,  lower  than  what 
is  expected. 

To  investigate  the  reason  for  this  discrepancy,  we  conducted  an  additional  experiment 
where  we  vary  the  normalized  rate  of  the  interferer  from  0  to  1,  to  determine  if  there  is  a 
better  rate  it  could  have  chosen.  Figure  8b  shows  the  received  rate  for  both  connections  as 
a  function  of  the  interferer’s  rate.  This  is  a  very  surprising  result,  because  the  rate  received 
by  the  lower  connection  is  non-linear,  contradicting  results  from  the  literature  [14],  This 
pattern  was  observed  with  other  pairs  of  nodes  in  our  experiments,  but  we  have  not  yet 
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Figure  7:  Wireless  mesh  topology 
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(a)  Performance  of  PI  model 


(b)  Non-linear  relationship  between  interferer  and  in- 
terferee  links 


Figure  8:  Experiments  with  partial  interference  model 
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Link  Sum  Utility 


Figure  9:  Sum  utility 


determined  the  exact  reasons  for  it.  This  is  a  major  area  for  future  research,  because  this 
may  require  significant  changes  to  our  underlying  model  of  the  network. 

Despite  this  finding,  the  partial  interference  model  does  perform  better  than  the  binary 
interference  model  on  this  topology,  in  terms  of  overall  link  utility.  Figure  9  shows  the 
link  utility  for  the  partial  interference  model  as  compared  to  binary  interference,  for  (a) 
when  interference  is  modeled  as  contention  and  (b)  when  interference  is  ignored.  Modeling 
interference  as  contention  is  far  too  conservative,  meaning  it  is  better  to  let  interfering  links 
corrupt  the  packets  of  other  links  rather  than  making  them  take  turns.  Ignoring  interference 
does  very  well,  due  to  the  non-linear  behavior  seen  in  Figure  8b.  However,  because  the 
interferer  is  not  rate  limited  at  all,  the  lower  connection  suffers  greatly  and  sometimes  is 
completely  starved.  Thus  even  though  the  partial  interference  model  is  not  exactly  correct, 
there  is  some  benefit  to  rate  limiting. 

42  First  Principles  Modeling  of  Wireless  Networks 

We  further  improve  the  accuracy  of  our  wireless  networks  model  by  focusing  on  the  con¬ 
straints  used  to  model  the  network,  as  this  is  the  critical  piece  in  the  optimization  approach. 
Here  we  use  a  first-principles  model,  meaning  that  we  start  with  the  most  basic  assumptions 
of  how  multi-hop  wireless  networks  with  CSMA  operate.  In  this  model,  perceived  times  that 
the  medium  is  occupied  are  represented  as  a  random  set.  Our  model  may  also  be  classified 
as  a  measurement  based  model,  as  it  takes  as  inputs  the  probabilities  of  finks  carrier  sens¬ 
ing  or  interfering  with  each  other.  Such  an  approach  is  more  realistic  than  physical  layer 
modeling,  such  as  the  various  signal  fading  models  with  SINR  thresholds,  and  has  been  used 
significantly  to  model  wireless  networks  for  various  purposes  [15,  8,  17,  16,  11,  14,  12], 

42.1  Model 

In  this  model,  we  make  the  following  elementary  assumptions: 

•  Discretization  of  time.  Time  is  divided  into  large  blocks  of  time  that  are  further  divided 
into  equally  sized  slots.  Dining  each  time  slot,  each  link  is  either  sends  or  doesn’t. 
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•  Uniform  random  selection.  For  each  time  block  T,  each  link  has  a  set  F  C  T  of 
available  time  slots  in  which  to  send,  and  a  set  X  C  F  when  it  does  send.  Each  t  e  F 
has  an  equal  probability  of  being  in  X. 

•  Negligible  indirect  scheduling.  Using  the  notation  from  the  above  assumption,  if  link  i 
carrier  senses  links  j  and  k  sending  during  Xj,  X^  C  T,  respectively,  then  dependencies 
of  Xj  U  Xk  on  the  sending  times  of  any  link  /  ^  i,j,  k  is  negligible. 

When  considering  the  effective  rate  of  links  j  and  k  as  perceived  by  link  i,  realistically 
there  may  be  some  other  link  /  (or  even  a  set  of  other  links)  that  cause  the  rates  of  j  and 
k  to  overlap  more  or  less  than  usual,  which  could  in  turn  affect  how  much  link  i  can  send. 
The  last  assumption  simply  states  that  these  effects  are  negligible.  We  recognize  that  it  can 
significantly  impact  the  accuracy  of  the  model.  This  concern  will  be  addressed  in  future 
research  with  empirical  testing. 


4.2.2  Problem  Formulation 

In  this  model,  the  sending  constraint  is  given  by 

Si  +  Si  <  1,  Vi  G  L, 

where 

Si=  (-t  )lphlfi(p)gi(p)h(p), 

pev(Li) 

fi(p)  =  II  chsT 
kip) 


9i(P )  = 


niepMY 

hip) = 1  -  *  yi  (-i)ip,hi  n 

p'ev(p)  jep1 

and  the  independence  is  given  by 

h(p)  J^J  (1  c%j  Cji  T  CijCji ). 
{i,j}eV2(p) 


The  receiving  constraint  is  given  by 


where 


and 


r  i  d{  ( 1  Ri  )si: 


Ri  =  (_1)|p|  1fi(p)h(p) 

P&V{Li) 

flip)  =  II  °ris:)- 
j&p 


(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 
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In  this  model,  V(p)  is  the  set  of  all  subsets  of  p  except  the  empty  set,  and  Vz(p)  is  the 
set  of  all  subsets  of  p  with  \p\  =  z. 

Note  that  (17)  is  an  approximation.  If  one  link  carrier  senses  another  link  completely 
(cij  =  1)  then  their  random  sets  do  not  intersect.  If  any  two  random  sets  in  p  do  not  intersect, 
then  the  intersection  of  p  is  empty,  which  means  that  h(p)  should  equal  zero.  Only  when  all 
random  sets  are  independent  should  it  equal  one. 

4.2.3  Numerical  Results 

Instances  of  this  problem  are  frequently  non-convex,  due  to  the  addition  and  subtraction 
of  several  rational  functions  in  Si,  and  due  to  similar  reasons  in  Ri.  We  have  developed  a 
branch  and  bound  solution  to  solve  the  problem  by  successively  dividing  the  hypercube  in 
which  the  feasible  set  resides  into  smaller  regions,  and  evaluating  lower  and  upper  bound 
functions  for  the  optimal  value  in  each  region.  The  bounds  on  each  region  allow  one  to 
conclude  that  some  regions  need  not  be  divided  further. 

We  use  this  solution  to  compare  the  performance  of  the  first-principles  model  against  the 
PI  and  BI  models  on  several  kinds  of  topologies.  To  do  so,  we  define  the  optimality  of  a 
controller  as 

O  =  P'/P*,  (21) 

where  P1  is  the  performance  (geometric  mean  of  receiving  rates)  of  the  controller  and  P* 
is  the  performance  (lower  bound)  reported  by  the  first-principles  NUM  problem.  Note  that 
O  is  simply  the  inverse  of  the  performance  ratio  R  used  earlier.  Thus  an  optimality  of  1 
equates  to  no  performance  loss  despite  the  inaccuracies  in  the  model  with  respect  to  partial 
carrier  sensing  or  partial  interference. 

For  the  smaller  topologies,  the  branch  and  bound  algorithm  converged  to  within  a  differ¬ 
ence  of  0.01  performance.  However,  for  some  of  the  larger  topologies,  the  algorithm  took  too 
long  to  converge.  We  therefore  will  also  report  a  certainty  measure  in  these  cases  according 
to  (21),  where  P'  is  the  branch  and  bound’s  lower  bound  score  and  P*  is  the  upper  bound 
score.  This  measure  tells  us  how  much  higher  the  true  optimal  score  might  be,  and  thus  how 
much  worse  the  optimality  of  the  controllers  might  be.  If  the  certainty  is  not  reported  for  a 
particular  topology,  then  it  was  very  close  to  1. 

The  binary  contention  graphs  for  the  partial  and  binary  interference  models  were  built 
based  on  the  independence  approximation  (17).  If  two  links  had  an  independence  greater 
than  0.5,  an  edge  was  drawn  between  them.  For  the  maximal  clique  model,  a  contention  edge 
also  needed  to  be  drawn  for  high  interference  values.  This  was  done  by  defining  a  function 
mimicking  the  independence,  but  using  a  instead  of  c.  An  edge  was  drawn  on  the  maximal 
clique  model’s  contention  graph  if  the  multiplication  of  the  two  “independence”  functions  ( a 
and  c)  was  greater  than  0.5. 

When  computing  results  over  a  range  of  interference  factor  values  a,  we  omitted  any 
topologies  such  that  there  existed  two  links  i  and  j  where 

o,j  ^  1  (‘ij  ■  (22) 

This  is  because  such  topologies  cannot  occur  due  to  the  relationship  between  a  and  c.  If  two 
links  carrier  sense  each  other  well,  their  sending  rates  cannot  overlap  much  and  therefore 
they  cannot  interfere  with  each  other  much. 
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(a)  Contention/interference  graph  (b)  Optimality  of  the  binary  inter-  (c)  Optimality  of  the  partial  inter- 

of  the  two-link  topologies  ference  controller  ference  controller 


Figure  10:  Optimality  of  the  controllers  for  the  two-link  topologies 


We  first  measured  optimality  on  the  simpliest  topologies  of  interest,  with  only  two  links. 
The  partial  contention/interference  graph  of  this  class  of  topologies  is  given  in  Figure  10a. 
In  this  figure  and  subsequent  figures,  we  let  solid  arrows  denote  contention  (where  the  arrow 
points  to  the  link  that  is  sensing)  and  dashed  arrows  denote  interference  (where  the  arrow 
points  to  the  link  being  interfered  with).  The  values  of  a12,  a2 1,  c12,  and  c21  were  permutated 
over  the  range  of  0  to  1,  with  granularity  of  0.2,  omitting  the  topologies  that  were  unrealistic. 
The  CDFs  of  the  optimality  for  both  models  appear  in  Figure  10b  and  Figure  10c.  This 
shows  that  the  binary  interference  controller  performs  above  0.9  optimality  in  most  cases, 
but  many  cases  drop  well  below  that.  The  partial  interference  controller,  on  the  other  hand, 
almost  always  performs  above  0.9  optimality. 

We  next  tested  optimality  over  various  three-link  topologies  by  choosing  the  two-link 
topology  with  the  worst  optimality  for  the  partial  interference  controller  and  adding  a  link  to 
it.  We  first  tested  over  all  interference  factors  a  both  from  and  to  the  third  link,  with  a  gran¬ 
ularity  of  0.2,  and  c  values  involving  link  3  set  to  zero.  The  partial  contention/interference 
graph  of  this  class  of  topologies  is  given  in  Figure  11a.  The  CDFs  of  the  optimality  for  both 
models  appear  in  Figure  lib  and  Figure  11c.  For  the  binary  interference  controller,  opti¬ 
mality  is  worst  at  0.56,  with  various  topologies  scoring  at  or  around  this  level.  Optimality 
is  close  to  1  largely  for  the  cases  where  there  is  high  interference  ( a  >  0.8)  for  most  or  all 
parameters.  For  the  partial  interference  controller,  95%  of  all  cases  are  at  0.85  optimality  or 
better. 

We  also  tested  over  all  contention  coefficients  c  both  from  and  to  the  third  link,  with  a 
granularity  of  0.2,  and  a  values  involving  link  3  set  to  zero.  The  partial  contention/interference 
graph  of  this  class  of  topologies  is  given  in  Figure  12a  and  the  optimality  CDFs  in  Figure  12b 
and  Figure  12c.  The  two  CDFs  are  identical  because  without  any  interference,  the  two  clas¬ 
sical  models  are  identical.  In  general,  these  controllers  lose  performance  when  there  is  a 
significant  amount  of  partial  carrier  sensing  in  the  network  (c  values  not  close  to  0  or  1). 
However,  we  have  seen  in  the  above  results  of  varying  the  a  parameters  that  introducing 
partial  interference  into  a  topology  plagued  with  partial  carrier  sensing  usually  results  in 
much  better  optimality  for  the  partial  interference  controller. 

Another  set  of  topologies  we  tested  over  was  with  /  interferers  on  one  link,  with  a  =  1 
and  c  =  0.5  between  the  interferers.  Figure  13a  shows  the  partial  contention/interference 
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(a)  Contention/interference  graph  (b)  Optimality  of  the  binary  inter-  (c)  Optimality  of  the  partial  inter- 

of  the  three-link  topologies  varied  ference  controller  ference  controller 

over  interference 


Figure  11:  Optimality  of  the  controllers  for  the  three-link  topologies  varied  over  interference 


(a)  Contention/interference  graph  (b)  Optimality  of  the  binary  inter-  (c)  Optimality  of  the  partial  inter- 

of  the  three-link  topologies  varied  ference  controller  ference  controller 

over  contention 


Figure  12:  Optimality  of  the  controllers  for  the  three-link  topologies  varied  over  contention 
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(a)  Contention/interference  graph 
of  the  topologies  with  partially  de¬ 
pendent  interferers 


(b)  Optimality  of  binary  and  par¬ 
tial  interference  models  for  par¬ 
tially  dependent  interferers 


Figure  13:  Optimality  of  the  controllers  for  partially  dependent  interferers 


graph  for  this  class  of  topologies.  Figure  13b  plots  the  optimality  of  the  classical  models 
across  values  of  I.  Only  when  7  =  2  does  the  partial  interference  controller  have  signif¬ 
icantly  low  optimality,  whereas  the  maximal  clique  controller  always  does  poorly.  These 
results  show  that  when  contention  between  interferers  is  only  partial,  the  partial  interference 
controller  performs  much  better  in  the  face  of  complete  interference  (a  =  1)  than  it  does 
when  interferers  are  completely  dependent. 

We  next  tested  optimality  on  a  chain  topology  of  5  nodes  arranged  in  a  line,  with  active 
links  sending  both  ways  between  each  node.  Since  many  of  the  links  are  connected  by  a  node, 
there  is  a  significant  amount  of  perfect  contention  (c  =  1),  typical  of  a  real  network.  There 
is  also  interference  among  nodes  that  are  farther  away  from  each  other.  On  this  topology, 
the  binary  interference  controller  achieved  an  optimality  of  0.861  and  the  partial  interference 
controller  achieved  an  optimality  of  0.978.  It  is  possible  that  the  true  optimalities  are  much 
worse,  since  the  certainty  of  the  first-principles  branch  and  bound  solution  is  only  0.709. 
However,  based  on  the  lengthy  runs  of  the  branch  and  bound  algorithm  we  believe  it  is 
unlikely  that  the  true  optimal  score  is  closer  to  the  upper  bound,  since  the  lower  bound 
never  moves  and  the  upper  bound  slowly  approaches  the  lower  bound.  The  branch  and 
bound  algorithm  ran  on  this  topology  for  approximately  16  hours  and  60,000  iterations. 

Finally,  we  tested  optimality  on  a  mesh  topolgy,  as  shown  Figure  14.  We  designed  this 
topology  to  be  similar  to  a  typical  flow  of  traffic  through  a  mesh  network  providing  Internet 
access  to  users,  and  is  based  on  the  positioning  of  nodes  on  one  floor  of  a  mesh  network 
deployed  in  our  department’s  building.  It  can  be  thought  of  as  a  multi-hop  Internet  access 
network  for  users,  where  nodes  X  and  Y  are  the  access  points  to  the  cloud.  Note  that  just 
about  any  subset  of  the  pairs  of  nodes  could  be  the  links  that  are  simultaneously  active,  and 
the  choice  for  this  particular  problem  is  somewhat  arbitrary.  The  two  parts  of  the  network 
do  contend  and  interfere  with  each  other. 

Like  the  chain  topology,  the  mesh  topology  has  a  good  mix  of  partial  and  binary  con¬ 
tention,  with  some  partial  interference.  The  maximal  clique  controller  obtained  0.924  op¬ 
timality,  the  partial  interference  controller  obtained  0.97  optimality,  and  the  branch  and 
bound  algorithm  returned  a  certainty  of  0.807  for  the  optimal  score,  after  running  for  several 
days  and  approximately  78,000  iterations.  As  a  side  note,  a  certainty  of  0.795  was  obtained 
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Figure  14:  Network  graph  of  the  mesh  topoogy 


after  only  a  few  minutes,  and  the  lower  bound  never  moved,  as  was  the  case  for  all  topologies 
tested. 

The  take-home  message  of  these  numerical  results  is  that  the  maximal  clique  controller 
usually  has  significant  performance  loss,  whereas  the  partial  interference  controller  seems  to 
only  have  significant  performance  loss  in  select  topologies  where  the  approximations  intro¬ 
duced  by  the  partial  interference  model  sometimes  cause  the  controller  to  either  significantly 
under-utilize  the  medium  or  significantly  starve  at  least  links.  Specifically,  we  have  seen  that 
when  a  network  is  plagued  with  a  lot  of  partial  carrier  sensing,  but  also  has  very  little  inter¬ 
ference,  the  partial  interference  model  treat  the  partial  carrier  sensing  as  full  carrier  sensing, 
and  thus  under-utilize  the  medium.  We  have  also  seen  that  when  a  network  has  several 
interferers  that  can  carrier  sense  each  other  (near)  perfectly,  and  when  that  interference  is 
very  strong  (nearly  1),  the  partial  interference  model  will  make  the  mistake  of  having  the 
clique  of  interferers  send  at  or  close  to  the  full  clique  capacity,  thus  starving  the  link  being 
interfered  with. 

Despite  these  significant  failures  in  specific  instances  for  the  partial  interference  model 
(optimalities  as  low  as  0.33),  the  results  on  the  larger  networks  of  the  chain  and  mesh 
topologies  indicate  that  perhaps  these  dangerous  motifs  either  do  not  occur  frequently,  or 
their  detriment  is  washed  out  by  the  good  performance  in  other  parts  of  the  network.  In 
both  topologies,  the  partial  interference  controller  scored  above  0.97.  Of  course,  the  results 
of  the  chain  and  mesh  topologies  cannot  prove  such  high  optimality  for  other  topologies  of 
similar  size — but  combined  with  the  near-exhaustive  search  over  the  small  topologies,  it  is  a 
strong  indication  that  the  partial  interference  model  is  good  enough  for  the  purposes  of  rate 
control. 

4.3  Robust  dynamical  network  structure  reconstruction 

One  of  the  fundamental  interests  in  systems  biology  is  the  discovery  of  the  specific  biochem¬ 
ical  mechanisms  that  explain  the  observed  behaviour  of  a  particular  biological  system.  In 
particular,  we  consider  the  problem  of  reconstructing  the  network  structure  from  input  and 
partially  measured  output  data  of  a  dynamical  system,  and  in  turn  uncovering  the  underlying 
mechanisms  responsible  for  the  observed  behaviour.  The  biological  network  reconstruction 
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problem  challenges  come  from  the  necessity  to  deal  with  noisy  and  partial  measurements  (in 
particular,  the  number  of  hidden/unobservable  nodes  and  their  position  in  the  network  is 
unknown)  taken  from  a  nonlinear  and  stochastic  dynamical  network. 

There  are  several  tools  in  the  literature  to  infer  causal  network  structures.  These  tools  are 
mainly  rooted  in  three  fields:  Bayesian  inference  [26,  40],  information  theory  (ARACNe  [21], 
[24],  [27])  and  ODE  methods  (inferelator  [22],  [19],  [20],  [28],  [37]).  Details  on  these  and 
other  methods  can  be  found  in  several  reviews  of  the  field  such  as  [18,  25,  30,  32,  36].  The 
vast  majority  of  network  reconstruction  methods  produce  estimates  of  network  structure 
regardless  of  the  informativity  of  the  underlying  data.  In  particular,  most  methods  pro¬ 
duce  estimates  of  network  structure  even  in  cases  with  data  from  only  a  few  experiments. 
Such  data  may  not  contain  enough  information  to  enable  the  accurate  reconstruction  of  the 
actual  network,  thus  the  obtained  network  estimates  can  be  arbitrarily  different  from  the 
true  network  structure  [25].  To  compensate  for  lack  of  information  in  data,  most  methods 
have  heuristics  that  try  to  “guess”  at  the  remaining  information,  either  by  specifying  prior 
distributions  or  by  appealing  to  a  priori  beliefs  about  the  nature  of  real  biological  networks, 
such  as  looking  for  the  sparsest  network.  Nevertheless,  these  heuristics  bias  the  results  and 
lead  to  incorrect  estimates  of  the  network  structure. 

In  contrast,  our  approach  has  been  to  identify  the  conditions  when  data  is  sufficiently 
informative  to  enable  accurate  network  reconstruction.  The  results  indicate  that  even  in 
an  ideal  situation,  when  the  underlying  network  is  linear  and  time-invariant  (LTI)  and  the 
measurements  are  noise-free,  network  reconstruction  is  impossible  without  additional  infor¬ 
mation  [29].  Surprisingly,  this  information  gap  is  not  due  to  a  lack  of  data,  or  a  deficiency 
in  the  number  of  experiments,  but  rather  it  occurs  because  system  states  are  only  partially 
observed;  the  information  gap  is  present  in  all  data  sets  except  those  that  satisfy  certain 
experimental  conditions.  Our  analysis  identified  a  particular  experimental  protocol  that  sat¬ 
isfies  these  necessary  conditions  to  ensure  that  data  will  be  sufficiently  informative  to  enable 
network  reconstruction.  This  protocol  suggests  that: 

1.  A  network  composed  of  p  measured  species  demands  p  experiments; 

2.  Each  experiment  requires  a  distinct  input  that  independently  controls  a  measured 
specie,  i.e.  experimental  input  i  must  affect  measured  specie  i  and  no  other  measured 
specie  except,  possibly,  indirectly  through  measured  specie  i. 

If  data  acquisition  experiments  are  not  performed  in  this  (or  an  equivalent)  way,  the  network 
cannot  be  reconstructed.  Moreover,  the  resulting  information  gap  is  catastrophic,  meaning 
that  any  internal  network  structure  explains  the  data  equally  well  (i.e.  fully  decoupled,  fully 
connected,  and  everything  in  between).  On  the  other  hand,  if  some  information  about  the 
network  is  available  a  priori,  as  is  usually  the  case,  then  these  conditions  can  be  relaxed  as 
explained  in  [29]. 

The  work  in  [29],  however,  did  not  take  into  account  the  realistic  scenario  that  typically 
systems  are  nonlinear  and  data  are  noisy.  This  section  extends  and  details  earlier  results  in 
[44]  by  developing  an  effective  method  to  reconstruct  networks  in  the  presence  of  noise  and 
nonlinearities,  assuming  that  the  conditions  for  network  reconstruction  presented  above  in 
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(1)  and  (2)  have  been  met.  Steady-state  (resp.  time-series)  data  can  be  used  to  reconstruct 
the  Boolean  (resp.  dynamical)  network  structure  of  the  system. 

The  section  is  organised  as  follows.  After  a  motivating  example  showing  that  input- 
output  data  alone  does  not  enable  network  reconstruction,  Section  7.3.1  reviews  dynamical 
structure  functions  and  gives  fundamental  results  concerning  their  usefulness  in  the  network 
reconstruction  problem.  Section  7.3.2  presents  the  main  results  of  the  section  regarding 
robust  network  reconstruction  from  input-output  data  subject  to  noise  and  nonlinearities. 
Finally,  we  conclude  this  section  with  biologically-inspired  examples  in  Section  7.3.3. 

Notation.  For  a  matrix  A  G  CMxJV,  AVJ  G  C  denotes  the  element  in  the  ith  row  and 
jth  column  while  Aj  G  CMxl  denotes  its  jth  column.  For  a  column  vector  a,  a[i]  denotes  its 
ith  element.  We  define  ejT  =  [0, . . . ,  0, 1  rth,  0, . . . ,  0]  G  Mlx7V.  /  denotes  the  identity  matrix. 
When  it  is  clear  from  the  context,  we  omit  the  explicit  dependence  of  transfer  functions  on 
the  Laplace  variable  s,  e.g.  we  write  G  instead  of  G(s). 

Motivating  example.  Consider  the  transfer  function 

i  r  i  i 


L  s+2  J 


G(s)  — 


s+ii 


obtained  from  data  (partial  observations)  using  system  identification  tools.  For  simplicity, 
assume  that  G(s)  accurately  represents  the  input-output  relation  of  the  original  system. 
This  transfer  function  is  consistent  with  two  state-space  realisations  x  =  Ax  +  Bu,  y  =  Cx 
given  by 
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Bi  —  B2  —  [0  0  1]T,  and  C\  =  C2  =  [. I  0]  G  M2x3  (i.e.,  the  third  state  is  hidden/non¬ 
observable).  Note  that  both  realisations  are  minimal  and  correspond  to  very  different  net¬ 
work  structures  as  seen  in  Figure  15.  This  demonstrates  that  even  in  the  idealised  setting 
(LTI  system,  non  noise  and  perfect  system  identification) ,  network  reconstruction  in  the 
presence  of  hidden/unobservable  states  is  not  possible  without  additional  information  about 
the  system. 


4.3.1  Dynamical  structure  functions  and  network  reconstruction 

In  [29]  we  introduced  the  notion  of  dynamical  structure  functions  and  showed  how  they  can 
be  used  to  obtain  necessary  and  sufficient  conditions  for  network  reconstruction.  For  the 
sake  of  clarity  and  completeness,  we  state  these  previously  obtained  results  here  without 
proofs.  We  refer  the  interested  reader  to  [29,  41]  for  the  corresponding  proofs. 

Consider  a  nonlinear  system  x  =  f(x,  u,  w\),  y  =  h(x,  w 2)  with  p  measured  states  y,  hid¬ 
den  states  z  (potentially  a  large  number  of  them),  m  inputs  u,  and  noises  wi,w2.  The  system 
is  linearised  around  an  equilibrium  point  (i.e.,  a  point  (x*,u*)  such  that  f(x*,u*,  0)  =  0), 
and  it  is  assumed  that  inputs  and  noises  do  not  move  the  states  too  far  from  the  equilibrium 
point  so  that  the  linearised  system  is  a  valid  approximation  of  the  original  nonlinear  system. 
The  linearised  system  can  be  written  as  x  =  Ax  +  Bu,  y  =  Cx,  where  x  =  x  — x*,u  =  u  —  u* 
and  y  =  h(x,  0)  —  h(x*,  0).  The  transfer  function  associated  with  this  linearised  system  is 
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Figure  15:  The  same  transfer  function  yields  two  minimal  realisations  with  very  different 
network  structures  (Left  vs.  Right).  Pink  nodes  are  measured,  while  blue  nodes  represent 
unmeasured  hidden  states;  the  top  diagram  on  either  side  reveals  the  complete  network  struc¬ 
ture  explicitly  showing  hidden  states,  while  the  lower  diagram  indicates  the  corresponding 
casual  structure  captured  by  the  dynamical  structure  function  (edges  associated  with  Q  are 
red,  while  those  associated  with  P  are  blue).  The  system  in  the  left  is  (kLi,  Bi,  C i)  in  7.3, 
and  the  system  on  the  right  is  (A2l  B2 ,  C2)-  Note  how  completely  different  the  two  network 
structures  are  (complete  decoupled  vs.  fully  connected)  even  though  either  realization  would 
be  an  equally  valid  description  if  all  one  knew  about  the  system  was  its  transfer  function, 
identified  from  input-output  data. 


given  by  G(s )  =  C(sl  —  A)  1B.  When  we  have  partial  observations,  i.e. ,  when  C  =  [I  0], 
we  partition  the  linearised  system  equation  as  follows 
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'  An 

7^12 
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where  x  =  [yT  zT]T  G  Mn,  is  the  full  state  vector,  y  G  is  a  partial  measurement  of  the 
state  (we  assume  p  >  1),  z  are  the  n  —  p  “hidden”  states,  and  u  G  Km  is  the  control  input. 
We  restrict  our  attention  to  situations  where  output  measurements  constitute  partial  state 
information,  i.e.,  p  <  n.  Taking  the  Laplace  transforms  of  the  signals  in  (23),  solving  for  Z, 
and  substituting  into  the  Laplace  transform  of  the  first  equation  of  yields  sY  =  WY  +  VU, 
where  W  =  An+A12  (si  —  A22)1  A2\  and  V  =  A12  (si  —  A22)~l  B2  +  Bi.  Now,  letting  D  be 
the  matrix  composed  of  the  diagonal  elements  of  W,  we  write  (si  —  D)Y  =  (W  —  D)  Y+VU . 
We  then  obtain  Y  =  QY  +  PU  where 

Q  =  (si  —  D)~l  (W  -  D)  and  P  =  (si  -  D)~l  V.  (24) 


Given  the  system  (23),  we  define  the  dynamical  structure  function  of  the  system  to  be  (Q,  P). 
If  all  the  measured  states  are  removed  from  the  system  except  for  Yt  and  Yj  then  the  transfer 
function  Qij  corresponds  to  the  exact  transfer  function  between  Yj  (considered  as  input)  and 
Yt  (considered  as  output).  The  same  holds  for  P  in  terms  of  Uj  and  Yt. 
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It  can  be  shown  that  G  =  (/  —  Q)~l  P  (see  [29]).  Based  on  this  latter  relation,  it  can  be 
seen  that  the  dynamical  structure  function  of  a  system  contains  more  information  than  the 
transfer  function,  and  less  information  than  the  state-space  representation  [29] .  We  can  then 
conclude  that,  with  no  other  information  about  the  system,  neither  dynamical  nor  Boolean 
reconstruction  is  possible.  Moreover,  for  any  internal  structure  Q  there  is  a  dynamical 
structure  function  ( Q,P )  that  is  consistent  with  G,  i.e.,  that  satisfies  G  =  (/  —  Q)~l  P.  In 
particular,  this  shows  that  the  use  of  criteria  such  as  sparsity  or  decoupledness  to  guide  our 
selection  of  a  proposal  network  structure  can  be  misleading.  If  one  were  to  optimise  for 
decoupledness,  for  example,  a  dynamical  structure  (0,  G)  could  and  would  always  be  found, 
regardless  of  the  true  underlying  structure.  Thus,  if  we  are  to  use  these  kinds  of  criteria, 
they  must  be  firmly  justified  a  priori. 

Proposition  1.  [29]  Given  a  pxm  transfer  function  G,  dynamical  structure  reconstruction 
is  possible  from  partial  structure  information  if  and  only  if  p  —  1  elements  in  each  column 
of  (■ Q ,  P)T  are  known  that  uniquely  specify  the  component  of  (Q,  P)  in  the  nullspace  of 

[GT  /]• 

The  importance  of  this  result  is  that  it  identifies  exactly  what  information  about  a  sys¬ 
tem’s  structure,  beyond  knowledge  of  its  transfer  function,  must  be  obtained  to  be  able  to 
recover  the  structure  without  appeal  to  a  priori  assumptions,  such  as  sparsity,  or  parsimony, 
etc.  This  enables  the  design  of  experiments  targeting  precisely  the  additional  information 
needed  for  reconstruction.  In  particular  when  p  =  m  and  G  is  full  rank,  we  observe  that 
imposing  that  P  is  diagonal,  i.e.,  that  each  input  controls  a  measured  state  independently, 
is  sufficient  for  reconstruction. 

Corollary  1.  [29]  If  m  =  p,  G  is  full  rank,  and  there  is  no  a  priori  information  about  the 
internal  structure  of  the  system,  Q,  then  the  dynamical  structure  can  be  reconstructed  if  each 
input  controls  a  measured  state  independently,  i.e.,  if,  without  loss  of  generality,  the  inputs 
can  be  numbered  such  that  P  is  diagonal. 

4.3.2  Robust  network  structure  reconstruction 

In  this  section,  we  consider  the  problem  of  robustly  reconstructing  dynamical  network  struc¬ 
tures.  Data  are  obtained  from  input-output  measurements  of  a  noisy  nonlinear  system. 
From  this  type  of  data  we  aim  to  find  the  internal  network  structure  Q  associated  with  the 
linearised  system  (23).  To  average  out  the  noise,  data-collection  experiments  are  repeated  N 
times.  For  simplicity  of  exposition,  we  assume  that  no  a  priori  information  on  the  internal 
network  structure  Q  is  available.  The  results  still  follow  if  some  a  priori  information  about 
Q  is  available,  and  such  information  can  typically  be  used  to  relax  the  experimental  protocol 
according  to  Proposition  1.  Hence,  data  are  collected  according  to  the  measurement  protocol 
described  in  the  introduction: 

(1)  the  number  of  distinct  data-collection  experiments  is  the  same  as  the  number  of  measured 
species.  This  in  particular  implies  that  u(t),y(t)  G  ; 

(2)  each  input  tp  controls  first  the  measured  state  yt  so  that  P  is  a  p  x  p  diagonal  matrix. 
In  the  following  two  subsections  (7.3.2  and  7.3.2),  we  propose  two  approaches  for  esti¬ 
mating  the  dynamical  structure  function  ( Q ,  P )  from  measured  input-output  data.  The  first 
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approach  is  indirect  and  involves  estimating  the  transfer  function  G  followed  by  computing 
(Q,  P)  from  G.  Since  some  information  is  lost  in  the  process  of  estimating  G,  we  consider  a 
second  approach  where  (Q,P)  is  directly  estimated  from  data  (without  estimating  first  G). 
Concerning  the  type  of  input-output  data  collected,  we  first  consider  time-series  input-output 
data  and  then  the  special  case  where  only  steady-state  data  are  available. 

Dynamical  network  reconstruction  from  identified  transfer  functions  This  section 
describes  a  method  to  obtain  the  dynamical  structure  function  from  a  transfer  function 
G.  This  transfer  function  was  identified  from  noisy  time-series  data  using  standard  system 
identification  tools  [33].  According  to  Corollary  1,  if  G  is  full  rank  there  is  a  unique  Q  and 
diagonal  P  satisfying  (/  —  Q)G  =  P.  Since  G  is  an  approximation  of  the  actual  system, 
Q  and  P  will  typically  be  mere  approximations  of  the  actual  dynamical  structure  function. 
Moreover,  due  to  noise  and  unmodelled  dynamics,  it  is  likely  that  Q  does  not  even  have  the 
correct  Boolean  structure.  Typically,  the  internal  structure  function  Q  obtained  from  such 
a  procedure  will  be  fully  connected,  i.e.,  all  non-diagonal  elements  of  Q  will  be  non-zero. 

The  main  idea  to  solve  the  network  reconstruction  problem  from  noisy  data  is  the  fol¬ 
lowing.  For  p  measured  states,  Q  has  p2  —  p  unknowns.  We  want  to  quantify  the  smallest 
distance  from  G  (or  directly  from  the  measured  data)  to  all  possible  Boolean  structures  (and 
there  are  2P  ~p  of  them).  Some  of  such  distances  will  be  large  revealing  that  the  correspond¬ 
ing  Boolean  structures  are  unlikely  to  be  the  correct  structures  while  other  will  be  small 
making  them  candidates  for  the  correct  structure. 

There  are  a  number  of  ways  to  model  input-output  data  with  noise  and  nonlinearities.  In 
order  to  obtain  a  convex  minimisation  problem,  we  consider  the  output  (could  also  be  input) 
feedback  uncertainty  model  [43].  In  this  framework,  the  “true”  system  is  given  by  (J+A)_1G, 
where  A  represents  unmodelled  dynamics,  including  nonlinearities,  and  noise.  Based  on  this 
choice  of  dynamic  uncertainty,  the  distance  from  data  to  a  particular  Boolean  structure  is 
chosen  to  be  ||A||,  in  some  norm,  such  that  Q  obtained  from  (/  +  A)_1G  =  (/  —  Q)~1P  has 
the  desired  Boolean  structure.  We  can  rewrite  the  above  equation  as  A  =  GP_1(J  —  Q)  —  I. 
Now,  let  X  =  P_1(J  —  Q).  Then  the  Boolean  structure  constraint  on  Q  can  be  reformulated 
on  X,  i.e.,  non-diagonal  zero  elements  in  X  correspond  to  those  in  Q  (since  Xl3  =  P^Qij 
for  i^j). 

We  can  order  all  Boolean  structures  from  1  to  2P“~P,  and  define  a  set  Xk  containing 
transfer  matrices  that  satisfy  the  following  conditions:  (i)  for  %  ^  j,  Xtj(s)  =  0  if  for  the 
considered  kth  Boolean  structure  Qij(s)  =  0;  all  other  Xij(s)  are  free  variables;  (ii)  when 
i  =  j,  Xu(s)  is  a  free  variable.  Hence,  the  distance  from  G  to  a  particular  Boolean  structure 
can  be  written  as  ak  =  inf'xeay  || GX  —  J||2,  which  is  a  convex  minimisation  problem  with 
a  careful  choice  of  a  norm.  Next,  we  show  that  this  problem  can  be  cast  as  a  least  squares 
optimisation  problem.  If  we  use  the  norm  defined  by  ||  A || 2  =  sum  of  all  ||Ajj|||,  where  ||  •  H2 
stands  as  the  ^2-norm  over  s  =  jui ,  then  using  the  projection  theorem  [39]  the  problem 
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reduces  to 


ak  =  inf  \\GX  -  If  =  inf  V  || GX,  -  e,f 

A  A  ‘ 
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=  Y  !vf  WAiYi~ei\\l 
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=  Y  WAWA^A^-etWl 

i 

where  Xt  is  the  ith  column  of  X  G  Xk,  Yi  is  a  column  vector  composed  of  the  free  (i.e., 
nonzero)  elements  of  Xtl  Ai  is  obtained  by  deleting  the  jth  columns  of  G  when  the  corre¬ 
sponding  elements  Xt[;j]  are  0  for  all  j,  and  (•)*  denotes  transpose  conjugate.  The  inhmum 
is  achieved  by  choosing  Y%  =  (A*Ai)~l A*ei,  and  A*A%  is  always  invertible  since  G  is  full  rank 

in  Corollary  1.  If  experiments  are  repeated  N  times,  yielding  a  transfer  function  G'  for  each 

r  iT 

experiment,  then  the  above  analysis  still  follows  simply  by  letting  G  =  GlT  ■  ■  ■  GnT 

Dynamical  network  reconstruction  directly  from  time-series  data  The  previous 
sections  used  a  two-step  approach  in  which  system  identification  was  first  used  to  estimate  a 
transfer  function  from  measured  input-output  data  and  then,  in  a  second  step,  the  identified 
transfer  function  was  used  to  obtain  a  dynamical  structure  function  representation  of  the 
system  which  is  optimal  in  terms  of  a  particular  metric.  This  section  proposes  a  method 
which  allows  identification  of  the  optimal  dynamical  structure  function  representation  di¬ 
rectly  from  the  measured  input  output  data.  The  advantage  of  this  direct  network  structure 
reconstruction  from  data  is  that  no  information  is  lost  during  the  initial  transfer  function 
identification  stage. 

Due  to  the  equivalence  between  dynamical  uncertainty  perturbations  [43],  we  are  free 
to  chose,  without  loss  of  generality,  the  type  of  uncertainty  perturbation  that  best  suits 
our  needs.  For  the  direct  method,  instead  of  a  feedback  uncertainty  as  was  considered  in 
the  previous  section,  the  uncertainty  perturbation  we  are  considering  here  is  the  additive 
dynamic  uncertainty  on  the  output,  i.e.,  Y  =  G^(U  +  A).  In  this  case,  we  think  about 
the  “distance”  in  terms  of  how  much  we  need  to  change  the  input  (data)  to  fit  a  particular 
Boolean  structure.  Since  G\  =  (/  —  Q)~lP  =  A'^1,  the  equality  Y  =  Ga{U  +  A)  can  be 
written  as 

A  =  XY  -  U, 

where  X  G  Xk,  for  some  particular  Boolean  network  k.  Recall  that  structural  constraints 
in  Q  can  be  imposed  directly  on  X  from  the  equality  X  =  P  l{I  —  Q).  We  can  therefore 
use  system  identification  tools  for  non-causal  autoregression  models  under  the  structural 
constraints  to  identify  X  (which  might  be  non-causal).  In  this  case,  the  distance  is  defined 
as  the  maximum  likelihood  of  the  estimation  problem. 

Penalising  connections  The  above  methodology  suffers  from  a  crucial  weakness:  there 
are  several  Boolean  structures  with  distances  smaller  or  equal  than  the  distance  to  the 
“true”  network.  Indeed,  the  extra  degrees  of  freedom  of  the  fully-connected  network  allow 
its  corresponding  distance  ak  to  be  the  smallest  of  all.  This  is  similar  to  the  noisy  data 
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over-fitting  problem  encountered  in  system  identification  where  the  higher  the  order  of  the 
transfer  function,  the  better  the  fit.  The  typical  approach  in  system  identification  is  to 
penalise  higher  dimensions  and  the  analogy  here  is  to  penalise  extra  network  connections. 

If  the  true  network  has  /  non-existent  connections  (/  off-diagonal  elements  in  Q  are  zero) 
then  there  are  2l  —  1  different  Boolean  networks  that  have  a  smaller  or  equal  distance  (due 
to  the  additional  degrees  of  freedom  provided  by  the  extra  connections).  When  noise  is 
present,  then  the  “true”  network  will  typically  have  an  optimal  distance  similar  to  these 
other  l  networks.  The  question  of  how  to  find  the  “true”  network  thus  arises.  With  repeated 
experiments,  small  enough  noise  (i.e.,  large  enough  signal-to-noise  ratio)  and  negligible  non- 
linearities,  the  optimal  distances  of  those  l  networks  are  comparable,  and  they  are  typically 
much  smaller  than  those  of  the  other  networks.  To  try  to  reveal  the  “true”  network,  one 
can  strike  a  compromise  between  network  complexity  (in  terms  of  number  of  connections) 
and  data  fitness  by  penalising  extra  connections.  There  are  several  ways  to  do  this.  Here, 
we  consider  one  of  the  classical  methods  known  as  Akaike’s  information  criterion  (AIC)  [31], 
or  some  of  its  variants  such  as  AICc  (which  is  AIC  with  a  second  order  correction  for  small 
sample  sizes),  and  the  Bayesian  information  criterion  (BIC)  [23]. 

The  AlC-type  approach  is  a  test  between  models  -  a  tool  for  model  selection.  Given  a 
data  set,  several  competing  models  may  be  ranked  according  to  their  AIC  value,  with  the 
one  having  the  lowest  AIC  being  the  best.  From  the  AIC  value  one  may  typically  infer  that 
the  best  models  are  in  a  tie  and  the  rest  are  far  worse,  but  it  would  be  arbitrary  to  assign  a 
threshold  above  which  a  given  model  is  rejected  [23].  The  AIC  value  for  a  particular  Boolean 
network  Bk  is  defined  as: 

AICk  =  2Lk-lnak,  (25) 

where  Lk  is  the  number  of  (non-zero)  connections  in  the  Boolean  network  Bk  and  ak  is  the 
optimal  distance  based  on  this  parameter  constraint. 

Although  finding  the  optimal  distance  in  the  second  term  of  eq.  (25)  can  be  done  effi- 

2 

ciently,  the  number  of  Boolean  networks  2P  ~p  grows  very  fast  with  the  number  of  measured 
states  p.  To  find  the  network  with  the  smallest  distance  it  is  thus  not  desirable  to  com¬ 
pute  the  optimal  distance  for  each  possible  Boolean  network.  Fortunately,  there  are  ways 
to  reduce  the  number  of  networks  that  need  to  be  considered.  As  we  saw  in  the  previous 
section  inixexk  [I dTAC  —  J||2  =  ]>A  inf^  ||  AiYj  —  e* 1 1 1  meaning  that  we  can  solve  each  optimi¬ 
sation  problem  separately.  Since  each  Fj  corresponds  to  p  —  1  unknowns  in  the  ith  row  of 
Q,  this  reduces  the  problem  to  solving  p2p_1  optimal  distances.  Finding  a  polynomial-time 
algorithm  to  compute  the  optimal  distance  through  this  method  is  a  subject  of  current  inves¬ 
tigation.  When  it  comes  to  the  steady-state  case,  [19]  proposed  a  polynomial-time  algorithm 
to  quickly  find  the  ranked  solutions  at  the  expense  of  solution  accuracy. 

Boolean  network  reconstruction  from  steady-state  data  So  far  we  have  assumed 
that  time-series  data  are  available.  Frequently,  however,  experimentation  costs  and  limited 
resources  only  permit  steady-state  measurements.  In  addition,  with  steady-state  measure¬ 
ments  it  is  typically  possible  to  perform  a  larger  number  of  experiments  within  the  same 
amount  of  time,  effort  and  cost.  As  shown  below,  most  of  the  connectivity  of  the  network 
together  with  the  associated  steady-state  gains  (and  the  associated  positive  or  negative  sign) 
can  still  be  reconstructed  from  steady-state  data.  However,  no  dynamical  information  will 
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be  obtainable.  In  other  words,  for  most  cases  we  can  still  recover  the  Boolean  network  from 
steady-state  data. 

Assume  that  after  some  time  of  maintaining  the  control  input  concentrations  at  a  constant 
value,  the  measured  outputs  y  have  converged  to  a  steady-state  value.  This  is  equivalent 
(if  the  system  is  stable  or  quasi-stable  [37])  to  assuming  that  we  can  obtain  G(0),  i.e.,  G(s) 
evaluated  at  s  =  0.  Now  the  relationship  (I  —  Q(s))G(s)  =  P(s )  evaluated  at  s  —  0  becomes 
(I—Q(0))G(0)  =  P( 0).  From  this  equation  and  the  knowledge  of  G(0),  all  of  the  results  given 
in  Sections  7.3.2  and  7.3.2  follow  provided  that  no  element  of  G(s)  has  a  system  zero  [43]  at 
0.  In  that  case,  a  nonzero  element  in  the  obtained  Boolean  network  indicates  the  existence  of 
a  causal  relationship  between  the  corresponding  pair  of  nodes  while  a  zero  element  indicates 
the  absence  of  such  relationship. 

4.3.3  Biologically-inspired  examples 

This  section  illustrates  with  two  examples  the  theoretical  results  presented  in  the  previous 
section.  The  corresponding  sets  of  ordinary  differential  equation  describing  the  dynamics  of 
the  considered  networks  are  used  to  generate  noisy  data,  which  are  then  fed  to  our  recon¬ 
struction  algorithm  in  order  to  assess  its  ability  to  recover  the  correct  network  structure. 

Single  feedback  loop  In  this  first  example,  we  consider  the  following  nonlinear  system: 


i)\  = 

V max 

Vl  +  K  +-3+“‘ 

Ivm  1  ^3 

(26) 

i)2  = 

— 2t/2  +  l-5zi  +  u2 

(27) 

y.3  = 

—  1.52/3  +  0.5^2  +  u3 

(28) 

Zi  = 

O.82/1  —  0.5^i 

(29) 

z2  = 

I.22/2  -  0.8^2 

(30) 

h  = 

Fly, 3  -  1.3z3 

(31) 

where  Vmax  =  0.5  and  Km  =  0.1.  Equation  (26)  includes  a  nonlinear  function  of  z3  known 
as  a  Hill  equation.  It  represents  a  negative  regulation  of  the  rate  of  reaction  of  ij\  by  Z3.  For 
simplicity,  all  other  terms  are  linear.  In  this  example,  p  =  3,  i.e.,  there  are  three  measured 
states  (jq,  y 2  and  y3 )  while  the  other  3  states  are  hidden  (zi,  z2  and  z3).  The  corresponding 
network  is  given  in  Figure  16(a). 

Three  experiments  were  performed.  In  each  experiment,  one  input  was  a  step  while 
the  others  were  set  to  zero  and  data  was  collected  for  each  of  the  measured  species.  The 
experiments  were  repeated  3  times  to  average  out  the  noise.  For  simplification,  in  this 
example,  only  steady-state  data  was  used.  Data  was  obtained  by  numerically  integrating 
the  differential  equations  in  (26)-(31)  and  adding  independent  Gaussian  noises.  The  ratios 
between  standard  deviations  and  means  of  the  steady-state  data  were  within  the  range 
[0.35, 1.15],  which  shows  that  noise  is  considerable. 

Since  the  true  network  has  3  elements  in  Q  equal  to  zero,  there  are  23  =  8  networks  with 
a  better  or  equal  optimal  cost.  Computing  the  corresponding  distances  and  AICc  values 
for  all  possible  Boolean  structures  between  the  three  measured  species,  we  observed  that 
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Figure  16:  (a)  Complete  network  with  all  the  states.  The  red  circles  represent  the  measured 
states  while  the  blue  circles  correspond  to  hidden  states,  (b)  Network  of  the  measured  states 
only. 


Figure  17:  (a)  Network  representing  the  dynamical  interaction  between  the  10  species  be¬ 
lieved  to  be  responsible  for  the  chemotactic  response  of  Rhodobacter  sphaeroides.  We  assume 
that  only  species  Y[f,  Y£  and  “motor”  are  measured  (circled  in  red),  (b)  Network  connecting 
the  measured  states  only. 


the  distance  decreases  by  an  order  of  magnitude  when  we  arrived  at  the  true  network.  In 
addition,  AIC,  BIC  and  in  particular  AICc  are  able  to  pick  the  correct  network. 

Chemotaxis  in  Rhodobacter  sphaeroides  This  section  considers  the  reconstruction 
of  the  biochemical  network  responsible  for  chemotaxis  in  Rhodobacter  sphaeroides.  The 
network  is  represented  in  Figure  17  (see  [35,  38]  for  a  detailed  explanation  of  this  model 
and  its  biological  interpretation).  It  involves  10  species  dynamically  interacting  through 
a  complex  set  of  interconnections.  To  illustrate  our  method,  consider  noisy  data  from  3 
species  only:  Yf,  Y£  and  the  “motor”  (circled  in  red  in  Figure  17(a)),  obtained  based  on 
simulations  of  the  nonlinear  ordinary  differential  equation  model  proposed  by  [35] .  We  follow 
our  prescribed  experimental  protocol  and,  for  simplification,  only  steady-state  data  are  used. 
Relatively  large  Gaussian  noise  was  added  to  the  collected  data  to  simulate  measurement 
noise  in  the  data  set. 

Based  on  the  complete  network  given  in  Figure  17(a),  the  correct  network  to  recover 
is  presented  in  Figure  17(b).  Computing  the  corresponding  distances  and  AICc  values  for 
all  the  26  =  64  possible  Boolean  networks,  we  observed  that  the  network  with  the  smallest 
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AICc  was  not  the  correct  network  in  Figure  17(b)  as  it  was  missing  the  Q 12  link.  A  closer 
look  at  the  noisy  steady-state  data  of  Kf  (from  a  step  input  in  ii2)  revealed  an  extremely 
large  ratio  between  its  standard  deviation  and  mean  value  (~  200),  showing  that  the  noise 
was  completely  overpowering  the  signal.  Indeed,  it  can  be  shown  that  Y(f  has  a  very  small 
influence  on  Yf  since  the  pathway  from  Y(f  to  Y J  includes  a  reversible  reaction  with  a  very 
small  rate  constant  (for  detail,  see  [42,  34]).  The  next  set  of  smallest  values  of  AICc  consists 
of  4  networks,  including  the  true  one.  If  necessary,  an  extra  experiment  can  be  performed 
to  further  discriminate  between  these  five  candidate  networks. 

4.3.4  Discussion 

This  section  proposes  a  new  network  reconstruction  method  in  the  presence  of  noise  and  non- 
linearities  based  on  dynamical  structure  functions.  The  key  idea  is  to  find  minimal  distances 
between  the  existent  data  and  the  data  required  to  obtain  particular  network  structures. 
The  method  was  illustrated  with  two  biologically-oriented  examples.  They  showed  that  even 
in  the  presence  of  nonlinearities  and  considerable  noise  network  reconstruction  was  possi¬ 
ble.  Eventually,  when  the  signal  to  noise  ratio  was  too  small,  reconstruction  was  no  longer 
possible,  but  that  is  true  irrespective  of  the  method  used. 

Obviously,  the  method  has  limitations  with  respect  to  nonlinearities.  With  stronger  non¬ 
linear  terms  eventually  the  method  fails.  For  example,  network  reconstruction  for  oscillatory 
systems  is  still  an  open  problem.  However,  when  applied  to  the  reconstruction  of  various 
equilibrium  point  models  given  in  the  literature,  we  observed  that  reconstruction  was  always 
possible  when  the  signal-to-noise  ratio  of  the  measured  data  was  not  too  small  (far  less  than 
!)■ 

A  final  note  regarding  the  application  of  this  methodology  to  real  data.  We  have  looked 
throughout  the  literature  for  real  data  and  none  of  the  available  data  that  we  found  satisfied 
the  conditions  necessary  for  accurate  network  reconstruction.  Some  problems  that  we  observe 
in  the  literature  include:  1)  many  publications  do  not  include  raw  data  (they  typically  only 
include  means  and  standard  deviations),  and  many  authors  indicate  that  they  no  longer  have 
their  data;  2)  some  microarray  data  do  not  include  repeats  and  others  were  obtained  using 
dual  channel  microarrays  that  only  give  ratios  between  channels,  making  it  impossible  to 
reliably  extract  gene  expression  intensities.  One  of  the  most  promising  papers  was  [25] ,  which 
followed  our  experimental  protocol,  and  their  raw  data  is  available.  We  found,  however,  that 
the  data  did  not  meet  the  conditions  necessary  for  network  reconstruction.  In  the  paper,  the 
authors  compared  different  network  reconstruction  methods  only  to  find  that  none  of  the 
methods  even  came  close  to  identifying  the  true  network.  Nevertheless,  because  the  authors 
compared  the  results  to  random  guessing,  they  report  that  “Reverse  engineering  based  on 
differential  equations  and  Bayesian  networks  correctly  inferred  regulatory  interactions  from 
experimental  data.”  We  disagree  with  their  conclusion,  and  point  out  that  because  the  data 
was  not  sufficiently  informative  in  the  first  place,  such  a  comparison  is  not  meaningful.  To 
better  understand  the  degree  of  the  lack  of  information  in  the  data,  we  considered  all  10 
subnetworks  consisting  of  3  nodes.  The  gap  in  distances  between  fully  decoupled  and  fully 
connected  networks  ranged  from  just  2%  to  a  maximum  of  70%.  This  shows  that  there  is  not 
enough  information  in  the  data  to  differentiate  between  Boolean  structures.  Note  that  this 
data  was  obtained  from  over-expression,  so  based  on  these  results,  we  hypothesise  that  over- 
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expressing  genes  saturates  the  translation  and  transcription  machinery,  making  linearisation 
a  poor  approximation  of  the  actual  system  dynamics.  Current  work  is  exploring  the  design 
of  experiments  on  known  systems  that  1)  satisfy  our  data  collection  protocol  to  ensure 
the  resulting  data  is  sufficiently  informative  for  network  reconstruction,  and  2)  facilitate  a 
comparison  of  various  methods  so  we  can  better  understand  how  different  techniques  perform 
in  situations  where  accurate  network  reconstruction  is,  in  fact,  possible. 

4.4  Minimal  realization  of  dynamical  structure  functions  and  its 
application  to  network  reconstruction 

Recently,  networks  have  received  an  increasing  amount  of  attention.  In  our  information- 
rich  world,  the  questions  of  network  reconstruction  and  network  analysis  become  crucial  for 
the  understanding  of  complex  systems  such  as  biological,  social,  or  economical  networks. 
In  particular,  the  analysis  of  molecular  networks  has  gained  significant  interest  due  to  the 
recent  explosion  of  publicly  available  high-throughput  biological  data.  In  this  context,  the 
question  of  identifying  and  analyzing  the  network  structure  at  the  origin  of  measured  data 
becomes  a  key  issue. 

In  some  occasions,  measured  data  is  given  in  the  form  of  input-output  time-series  that 
describes  the  effect  of  inputs  on  outputs  (measured  states)  of  a  network.  When  data  is 
generated  by  a  linear  system,  a  matrix  transfer  function  describing  the  dynamic  input- 
output  behavior  is  generally  obtained  using  system  identification  [33].  If  the  original  state- 
space  model  is  available  or  deducible,  then  the  associated  network  structure  can  be  readily 
obtained  from  it.  However,  a  transfer  function  cannot,  in  general,  recover,  or  realize,  the 
original  state-space  model  since  the  realization  problem  does  not  typically  have  a  unique 
solution,  i.e.,  different  state-space  realizations  can  generate  the  same  input-output  behavior. 
Since  each  of  these  realizations  may  suggest  entirely  different  network  structures,  it  is  in 
general  impossible  to  identify  network  structures  from  transfer  functions  alone.  Therefore, 
more  information,  beyond  input-output  data  used  to  identify  a  transfer  function,  is  needed 
to  prefer  one  state-space  realization  over  another  as  a  description  of  a  particular  system  [45]. 

Another  difficulty  in  the  network  reconstruction  problem  comes  from  the  fact  that  the 
realization  problem  becomes  ill  posed  when  some  of  the  states  are  unobservable  or  “hidden” 
(this  even  happens  with  just  one  hidden  state  [43,  pp.  78]).  As  a  result,  failure  to  explicitly 
acknowledge  the  presence  of  hidden  states  and  the  resulting  ambiguity  in  network  structures 
can  lead  to  a  deceptive  and  erroneous  process  for  network  structure  discovery.  Consequently, 
determining  from  measured  data  the  presence  or  absence  of  a  causal  relationship  between 
two  variables  in  a  network  is  a  challenging  question. 

Motivated  by  this,  we  are  focusing  on  the  effect  of  hidden  states  in  the  network  that  we 
are  aiming  to  reconstruct.  A  new  representation  for  LTI  systems,  called  dynamical  structure 
functions  was  introduced  in  [29].  Dynamical  structure  functions  capture  information  at  an 
intermediate  level  between  transfer  function  and  state  space  representation  (see  Figure  18). 
Specifically,  dynamical  structure  functions  not  only  encode  structural  information  at  the 
measurement  level,  but  also  contain  some  information  about  hidden  states.  Based  on  the 
theoretical  results  presented  in  [29],  we  proposed  some  guidelines  for  the  design  of  an  ex¬ 
perimental  data-acquisition  protocol  which  allows  the  collection  of  data  containing  sufficient 
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Figure  18:  Mathematical  structure  of  the  network  reconstruction  problem  using  dynamical 
structure  functions.  Red  arrows  mean  “uniquely  determine” ,  blue  arrows  indicate  our  work. 


information  for  the  network  structure  reconstruction  problem  to  become  solvable.  In  par¬ 
ticular,  we  have  shown  that  if  nothing  is  known  about  the  network,  then  the  data-collection 
experiments  must  be  performed  as  follows: 

(A.l)  for  a  network  composed  of  p  measured  species,  the  same  number  of  experiments  p  must 
be  performed; 

(A. 2)  each  experiment  must  independently  control  a  measured  species,  i.e. ,  control  input  i 
must  first  affect  measured  species  i. 

If  the  experiments  are  not  performed  in  this  way  the  network  cannot  be  reconstructed,  and 
any  network  structure  fits  the  data  equally  well  (e.g.  a  fully  decoupled  network  or  a  fully 
connected  network).  If  biologists  have  already  some  information  about  the  network,  as  it  is 
usually  the  case,  then  these  conditions  can  be  relaxed  as  explained  in  [29]. 

Using  dynamical  structure  functions  as  a  mean  to  solve  the  network  reconstruction  prob¬ 
lem,  the  following  aspects  need  to  be  considered  (see  Figure  18): 

First  (see  (A)  in  Figure  18),  the  properties  of  a  dynamical  structure  function  and  its 
relationship  with  the  transfer  function  associated  with  the  same  system  need  to  be  precisely 
established  (this  was  done  in  [29]). 

Second  (see  (B)  in  Figure  18),  an  efficient  method  is  developed  to  reconstruct  networks 
in  the  presence  of  noise  and  nonlinearities  (this  was  done  in  [46]).  This  method  relies  on 
the  assumption  that  the  conditions  for  network  reconstruction  presented  above  in  (A.l)  and 
(A. 2)  have  been  met.  In  our  approach,  we  use  the  same  information  as  traditional  system 
identification  methods,  i.e.,  input-output  data.  However,  with  our  method,  steady-state 
(resp.  time-series  data)  can  be  used  to  reconstruct  the  Boolean  (resp.  dynamical  network) 
structure  of  the  system  (see  [46]  for  more  details). 

Third  (see  (C)  in  Figure  18),  once  the  dynamical  structure  function  is  obtained,  as 
a  main  result  of  this  section,  an  algorithm  for  constructing  a  minimal  order  state-space 
representation  consistent  with  such  function  is  developed.  In  an  application,  this  provides 
a  way  to  estimate  the  complexity  of  the  system  by  determining  the  minimal  number  of 
hidden  states  in  the  system.  For  example,  in  the  context  of  biology  it  helps  understand  the 
number  of  unmeasured  molecules  in  a  particular  pathway:  a  low  number  means  that  most 
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molecules  in  that  pathway  have  been  identified  and  measured,  showing  a  good  understanding 
of  the  system;  while  a  large  number  shows  that  there  are  still  many  unmeasured  variables, 
suggesting  that  new  experiments  should  be  carried  out  to  better  characterize  that  pathway. 

The  outline  of  the  section  is  as  follows.  Section  7.4.1  reviews  the  definition  of  dynamical 
structure  functions  and  their  properties.  The  main  result  can  be  found  in  Section  7.4.2 
where  we  propose  a  minimal  order  realization  algorithm  based  on  state-space  realizations 
and  pole-zero  analysis.  Simulation  and  discussion  are  addressed  in  Section  7.4.3.  Finally 
conclusions  are  presented  in  Section  7.4.4. 


4-4.1  System  Model 

Consider  a  linear  system  (it  can  also  be  a  linearization  of  some  original  nonlinear  system) 
x  =  Ax  +  Bu ,  y  =  Cx.  The  transfer  function  associated  with  this  system  is  given  by 
G(s)  =  C(sl  —  A)^1B.  Typically,  we  can  use  standard  system  identification  tools  [33]  to 
identify  a  transfer  function  G(s)  from  input-output  data. 

Like  system  realization,  network  reconstruction  also  begins  with  the  identification  of  a 
transfer  function,  but  it  additionally  attempts  to  determine  the  network  structure  between 
measured  states  without  imposing  any  additional  structure  on  the  hidden  states.  As  we 
have  shown  in  [29],  this  requires  a  new  representation  of  linear  time- invariant  systems:  the 
dynamical  structure  function  (defined  later).  An  algorithm  allowing  the  dynamical  structure 
function  to  be  obtained  from  input-output  data  is  proposed  in  [46].  This  section  assumes 
that  the  dynamical  structure  function  has  already  been  obtained  from  data,  and  will  focus 
on  finding  one  of  its  minimal  state-space  realizations. 

The  dynamical  structure  function  is  obtained  as  follows:  First,  we  transform  [A,  B,  C\  to 
[A°,  B°,  [Jp  0]]  without  changing  G(s),  where  p  =  rank(C').  The  linear  system  dynamics 
then  writes 

y 

z 

y 

where  x  =  ( y ,  z )  G  Mn°  is  the  full  state  vector,  y  EW  is  a  partial  measurement  of  the  state,  z 
are  the  n°  —p  “hidden”  states,  and  u  G  Mm  is  the  control  input.  In  this  work  we  restrict  our 
attention  to  situations  where  output  measurements  constitute  partial  state  information,  i.e., 
p  <  n°.  We  consider  only  systems  with  full  rank  transfer  functions  that  do  not  have  entire 
rows  or  columns  of  zeros,  since  such  “disconnected”  systems  are  somewhat  pathological  and 
only  serve  to  complicate  the  exposition  without  fundamentally  altering  our  conclusions. 

Taking  the  Laplace  transforms  of  the  signals  in  (32)  yields 
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where  Y,  Z,  and  U  are  the  Laplace  transforms  of  y,  z,  and  u,  respectively, 
gives 


Z  =  (sl-  A^)-1  A21Y  +  (si  -  A^)"1  B°2U 


Solving  for  Z 
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Substituting  this  last  expression  of  Z  into  (33)  then  yields 

sY  =  W°Y  +  V°U  (34) 

where  W°  =  A°n  +  A°l2  {si  —  A^y1  A21  and  V°  =  B°  +  A°l2  {si  —  A22)~L  B2. 

Now,  let  R°  be  a  diagonal  matrix  formed  of  the  diagonal  terms  of  W°  on  its  diagonal, 
i.e.,  R°  =  diag{H/0}  =  diag(hh1°1,  W22, ...,  W°p).  Subtracting  R°Y  from  both  sides  of  (34),  we 
obtain: 

{si  -  R°)  Y  =  {W°  -  R°)  Y  +  V°U 

Note  that  W°  —  R°  is  a  matrix  with  zeros  on  its  diagonal.  We  thus  have: 

Y  =  QY  +  PU  (35) 

where 

Q  —  {si  —  R T1  (r  -  R°) 

and 

p  =  {si  -  R°yl  v° 

Note  that  Q  is  zero  on  the  diagonal. 

Definition  1.  Given  the  system  (32),  we  define  the  dynamical  structure  function  of  the 
system  to  be  [Q ,  P] . 

Note  that,  in  general,  Q(s)  and  P{s)  carry  a  lot  more  information  than  G{s).  This  can 
be  seen  from  the  equality  G(s)  =  {I  —  Q{s))^1P{s)  (see  [29]  for  details).  However,  Q{s)  and 
P{s)  carry  less  information  than  the  state-space  model  (32)  (see  [46]). 

Definition  2.  A  dynamical  structure  function,  [Q,P],  is  said  to  be  consistent  with  a  par¬ 
ticular  transfer  function,  G,  if  there  exists  a  realization  of  G,  of  some  order,  and  of  the 
form  (32),  such  that  [ Q,P ]  are  specified  by  (36)  and  (37).  Likewise,  a  realization  is  consis¬ 
tent  with  [Q,P]  if  that  realization  gives  [Q,  P]  from  (36)  and  (37). 

Definition  3.  We  say  that  a  realization  is  G  minimal  if  this  realization  corresponds  to  a 
minimal  realization  of  G.  We  say  that  a  realization  is  [Q,P]  minimal  if  this  realization 
is  consistent  with  [Q,  P]  and  its  order  is  smaller  than  or  equal  to  that  of  all  realizations 
consistent  with  [Q ,  P] . 

The  underlying  principle  to  find  a  [Q,  P]  minimal  realization  is  to  search  for  a  realization 
with  the  minimal  number  of  hidden  states.  Such  a  realization  is  characterized  by  the  minimal 
number  of  pole- zero  cancellations  in  the  transfer  functions  Q  and  P. 

Proposition  2.  Given  a  dynamical  system  (32)  and  the  associated  dynamical  structure  func¬ 
tions  [Q,P]  with  R°  constructed  as  explained  above  (see  (32)-(37)j,  the  following  conditions 
must  hold 


(36) 

(37) 


A 


o 

11 


diag{A°u} 

diag{A°u} 

B° 


lim  R°(s); 

(38) 

s— >oo 

lim  sQ{s); 

(39) 

s— >•  oo 

lim  sP{s). 

s — S'-OO 

(40) 
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Proof.  Eq.  (38)  is  directly  obtained  from  the  definition  of  R°(s ): 

lim  R°(s)  =  lim  diag{hE°(s)} 

S— >•  OO  S— >•  OO 

=  diag{  lim  W°(s)}  =  diag-fA^} 

s — ^OO 


Since  the  proofs  for  eq.  (39)  and  (40)  are  very  similar,  we  focus  on  eq.  (39)  only.  Using  the 
fact  that  for  any  square  matrix  M,  (/  —  M)~l  =  YZo^1,  we  obtain,  from  the  definition  of 
Q  given  in  (36),  Q(s)  =  YZis~iR°  -  R°)  and  W°  =  A°u  +  YZ A°viKtX  Kx- 

Hence,  Q(s)  =  (A°±  —  R°(s))s~ 1  +  r(s),  in  which  r(s)  is  a  matrix  polynomial,  whose  largest 
degree  is  —2.  Finally,  multiplying  by  s  on  both  sides  and  taking  the  limit  as  s  goes  to  oo 
results  in  eq.  (39).  A  similar  argument  can  be  used  to  prove  eq.  (40).  □ 


We  give  an  illustrative  example  here. 


Example  1.  Consider  a  system  with  the  structure  depicted  in  Fig.  19.  A  linear  system’s 
representation  is 


an 

0 

®13 

0 

0  ' 

hi 

0  ' 

0 

®22 

0 

®24 

0 

0 

b-22 

0 

«32 

«33 

0 

«35 

X  + 

0 

0 

041 

0 

0 

®44 

0 

0 

0 

_  0 

0.52 

0 

0 

^55_ 

0 

0 

y  =  [I3  0]  x 


where  I3  is  the  3x3  identity  matrix.  Following  the  definitions  in  (36)  and  (37),  we  can 
write  down  the  corresponding  dynamical  structure  function  [Q,  P]  with 


Q 


p 


(  0 

^24^41 

(s  —  CL22){s  —  <244) 

V  0 


0 

0 

Q.35a52+a.32(s— a.55) 
(s-a33)(s-a55) 


a  13  \ 
s-an 

0 

0  J 


To  illustrate  Proposition  2,  we  have 


/0  0 

lim  sQ(s)  —  (  0  0 

s->°o  \  n 

\0  a32 


lim  sP(s)  =  s 

s— >00 


^>11 

s— an 


0 

0 


Generally,  there  exist  many  realizations  consistent  with  [Q,P],  In  the  following  section, 
we  focus  on  finding  a  [Q,P]  minimal  realization  (A,  B,  \j  0]),  i.e.,  a  realization  which  is 
consistent  with  [Q,  P]  and  which  has  minimal  order,  i.e.,  with  the  dimension  of  A  minimal 
(and  hence  the  lowest  possible  complexity). 
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(a)  (b) 


Figure  19:  (a)  An  example  system  with  two  inputs,  three  measured  states  (states  1,  2,  and  3) 
and  two  hidden  states  (states  4  and  5).  (b)  The  corresponding  dynamical  network  structure. 


4.4.2  Algorithm  to  find  a  [Q,  P }  minimal  realization 

From  a  dynamical  structure  function  [Q,  P]  we  cannot  reconstruct  [W°,  V°]  since  there  is  no 
information  regarding  the  diagonal  transfer  function  matrix  R°.  Given  [Q,  P]  and  a  diagonal 
proper  transfer  function  matrix  R ,  a  minimal  realization  of  [W  V]  =  [(si  —  R)Q  +  R  (si  — 
R)P]  can  be  obtained  as  follows: 

[W  V]  =  [An  B1]  +  A12(sI-A22)~1[A21  B2 ]  (41) 

The  idea  is  to  start  with  an  arbitrarily  chosen  R ,  and  then  use  a  state-space  realization 
approach  to  find  a  R*  which  minimizes  the  order  of  a  minimal  realization  of  [ W  V]. 

Lemma  1.  Suppose  W,  A22  and  A  are  defined  as  in  eq.  (41),  then  V  and  G  share  the  same 
zeros. 


Proof.  Since  si  —  W  is  the  Schur  complement  of  si  —  A22  in  si  —  A,  then 

det (si  —  A) 


det (si  —  W)  = 


det (si  -  A22) 


(42) 


Recall  that  V  =  B1  +  Ai2(sl  —  A22)  1B2.  Since  (si  —  W)G  =  V,  we  thus  have  that  V  and 
G  share  the  same  zeros  [49,  page  153].  □ 

Given  a  dynamical  structure  function  [Q,  P ],  a  random  choice  of  a  proper  diagonal  trans¬ 
fer  function  matrix  R  is  likely  to  result  in  additional  zeros  in  V  =  (sI  —  R)P.  From  Lemma  1, 
this  will  lead  to  additional  zeros  in  G  which  are  associated  to  uncontrollable  eigenvalues  of 
the  considered  realization  [47,  Section  4]  and  of  course  does  not  lead  to  a  minimal  realization 
in  eq.  (41).  At  this  stage  the  following  question  arises:  how  can  we  find  a  proper  diagonal 
transfer  function  matrix  R*  such  that  a  minimal  realization  of  [W  V ]  is  a  [Q ,  P]  minimal 
realization,  i.e., 


R*  =  argmin^  deg  {(si  —  R)s  1[sQ  sP]  +  [R  0]}  ,  (43) 

where  deg  is  the  McMillan  degree  [43].  Note  that,  since  there  are  many  choices  for  R*  that 
minimize  the  order  of  minimal  realizations  of  [ W  V],  a  chosen  R*  may  be  different  from  R° . 
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Assume  that  all  elements  in  [Q  P]  only  have  simple  poles.  This  assumption  can  be 
relaxed  but  we  adopt  it  here  for  simplicity.  Also  assume  that  [Q  P]  does  not  possess  any 
poles  at  0  (otherwise  we  can  change  eq.  (43)  to  (si  —  R)(s  —  a)-1[(s  —  a)Q  (s  —  a)P]  +  [F  0], 
where  a  G  R  is  not  a  pole  of  [Q  P]). 

Proposition  3.  Assume  [I  —  Q  P]  only  has  simple  poles  and  does  not  have  any  zeros.  A 
minimal  order  realization  of  [W  V]  in  (41)  can  be  achieved  using  a  constant  diagonal  matrix 

R*. 

Proof.  Assume  R*  has  at  least  one  term  on  the  diagonal  with  the  degree  of  numerator 
greater  or  equal  to  1,  e.g.,  suppose  the  ith  term  in  (si  —  R*)s _1  =  with  any  6  G  K 

and  deg(e.i(s ))  =  deg(<j>i(s ))  >  1,  where  deg(-)  returns  the  degree  of  a  polynomial.  Hence, 
the  multiplication  (si  —  i?*)s_1[sQ  sP]  will  introduce  deg(f>i(s ))  new  poles  and,  due  the 
assumption  of  simple  poles,  can  at  most  eliminate  deg(ei(s ))  =  deg(<fi(s ))  poles.  As  a 
consequence,  we  can  change  the  ith  term  to  ^4^  without  increasing  the  order.  Doing  this 
along  all  the  elements  of  R*  proves  the  result.  □ 

If  R*  is  a  constant  matrix,  the  term  [R*  0]  in  eq.  (43)  is  also  a  constant  matrix.  Therefore, 
the  order  of  a  minimal  realization  is  only  determined  by  (si  —  F*)^1)^  sP]  =  N[sQ  sP ]. 
Thus,  finding  the  “optimal”  R*  which  leads  to  the  minimal  order  in  eq.  (43)  is  equivalent 
to  finding  a  diagonal  proper  transfer  matrix  N  (N  with  corresponding  minimal  realization 
(A2,  B-2,  C2, 1)  is  restricted  to  the  set  of  matrices  of  the  form  (si  —  F*)s_1  with  a  constant 
R*  from  Proposition  3)  such  that  W[sQ  sP]  has  as  few  poles  as  possible.  Based  on  this 
idea,  the  following  algorithm  is  proposed: 

Step  1  :  Find  a  Gilbert ’s  realization  of  the  dynamical  structure  function. 

First,  using  the  results  in  [29,  Lemma  1],  we  find  a  minimal  realization  (Ai,  Bi,C\,  Di)  of 
[.sQ  sP] .  When  [sQ  sP]  has  /  simple  poles,  using  Gilbert’s  realization  [48]  gives 

1  TV. 

[sQ  sP }  =  +  lim  [s<2  sP], 

Z J  S  —  An  S— >00 

i=  1 

where  Kt  =  linis^At(s  —  A:)[sQ  sP]  and  has  rank  1  since  we  are  assuming  that  [sQ  sP]  has 
simple  poles. 

Consider  a  matrix  decomposition  of  K,  in  the  following  form: 

Ki  =  EiFi,  Vi, 

where  E%  e  and  F,:  =  (Ej Ef)^1  Ef  Ki.  Then  Ai  =  diag{A*}  e  M.lxl, 

Bx  =  [if  F2t  ...  P/r]T,  C\  =[E1  E2  ...  E,]  and  D1  =  linWoo[SQ  sP}. 

Step  2  :  Find  the  maximal  number  of  cancelled  poles. 

We  define  as  a  largest  subset  of  {El}  •  •  •  ,  such  that  all  the  elements  in  are  mutually 
orthogonal.  We  also  define  as  the  cardinality  of  <L.  Computationally,  (f  can  be  obtained 
using  the  algorithm  presented  in  the  Appendix.  We  claim  that  <f>  is  equal  to  the  maximum 
number  of  poles  we  can  eliminate  (the  proof  is  in  the  Appendix).  Therefore,  the  minimal 
order  of  [W  V }  is 

l  —  (j). 


APPROVED  FOR  PUBLIC  RELEASE:  DISTRIBUTION  UNLIMITED 

36 


As  a  consequence,  the  order  of  the  minimal  reconstruction  is  the  dimension  of  An  (constant 
p )  plus  the  minimal  dimension  of  A2 2  (obtained  above):  p  +  l  —  (j). 

Step  3  :  Construct  R*  to  obtain  the  minimal  reconstruction. 

Once  we  have  <h,  we  know  that  N(Xi)[j,j]  =  0  implies  R*[j,j]  =  A*.  Consequently,  each 
element  in  the  set  <f>  will  determine  at  least  one  element  in  R*.  This  last  fact  can  be  used  to 
construct  R*  element  by  element.  Once  R*  is  found,  we  can  obtain  A,  B  using  eq.  (41). 


4.4.3  Illustrative  example 

Example  2.  Consider  a  dynamical  structure  function  [Q ,  P] : 

\Q  I  p]  = 

We  first  compute  the  McMillan  degree  of  the  corresponding  transfer  function:  deg{G}  = 
deg  {(I  —  Q)_1P)}  =  4,  meaning  that  a  4th  order  state- space  model  is  enough  to  realize  the 
transfer  function.  It  is  interesting  to  see  what  is  the  minimal  order  realization  consistent 
with  the  dynamical  structure  function.  The  different  steps  of  the  algorithm  proposed  in  the 
previous  section  successively  yield  the  following: 

Step  1:  A  minimal  Gilbert  realization  of  s[Q,P]  is 


"  0 

1 

1 

1  - 

s+2 

s+3 

s+4 

1 

0 

1  1 

1 

s+1 

s+3 

s+4 

1 

1 

0 

1 

-s+1 

s+2 

s+4. 

M  =  diag{- 1,  -2,  -3,  -4},  P,  =  diag{ 2,  2,  2, 4}, 


0  -1  -1.5  -1 

0  111 

CA  = 

-0.5  0  -1.5  -1 

,  Dx  = 

10  11 

-0.5  -1  0  -1 

110  1 

Step  2:  By  definition,  Et  =  C\vt  where  Vi  G  M4  has  1  in  its  ith  position  and  zero  otherwise. 
Thus, 


{Ei,  •  •  •  ,  Efi\ 


Furthermore,  is  1  and  the  order  of  a  minimal  realization  of  the  given  dynamical  structure 
function  isp  +  l  —  fi  =  3  +  4  —  1  =  6.  Hence,  the  system  must  contain  at  least  3  hidden  states. 

Step  3:  R*  can  be  chosen  as  diag{a,  — 1,-1},  diag{— 2,  a,  —  2},  diag{  —3.  —3.  a},  or 
diag{—  4,  —4,  —4}  for  any  a  G  M. 

The  reconstructed  networks  are  represented  in  Fig.  20.  There  are  three  measured  (red) 
nodes,  labeled  1,2,3  and  by  the  analysis  above,  there  are  at  least  three  hidden  nodes  such 
that  the  corresponding  realization  is  consistent  with  the  dynamical  structure  function.  The 
red  connections  between  measured  nodes  are  the  same  for  all  candidate  networks  which  is  in 
accordance  with  Proposition  2.  Dashed  lines  correspond  to  the  connections  between  hidden 
and  measured  nodes. 

From  a  biological  perspective,  this  indicates  that  there  are  at  least  3  unmeasured  species 
interacting  with  the  measured  species.  Of  course,  the  “true”  biological  system  might  be  even 
more  complicated,  i.e.,  it  might  have  more  than  6  species.  Yet,  when  more  states  are  mea¬ 
sured,  the  dynamical  structure  functions  can  be  easily  updated  and  a  new  search  for  a  minimal 
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Figure  20:  Topologies  corresponding  to  the  four  [Q,  P]  minimal  realizations.  The  measured 
nodes  are  colored  red,  while  the  hidden  ones  blue.  Red  connections  between  measured  nodes 
are  the  same  for  all  the  networks  due  to  Proposition  2.  Each  node  has  a  self-loop  but  we 
omit  it  for  simplicity. 


realization  of  the  updated  system  can  be  performed  to  reveal  the  corresponding  minimal  num¬ 
ber  of  hidden  states. 

4.4.4  Summary 

In  this  section,  we  have  presented  a  method  for  obtaining  a  minimal  order  realization  consis¬ 
tent  with  a  given  dynamical  structure  function.  We  show  that  the  minimal  order  realization 
of  a  given  dynamical  structure  function  can  be  achieved  by  choosing  a  constant  diagonal  ma¬ 
trix  R*.  This  provides  a  way  to  estimate  the  complexity  of  the  system  by  determining  the 
minimal  number  of  hidden  states  that  needs  to  be  considered  in  the  reconstructed  network. 
For  example,  in  the  context  of  reconstruction  of  biological  networks  from  data,  it  helps  to 
understand  the  minimal  number  of  unmeasured  molecules  in  a  particular  pathway. 

4.5  The  Meaning  of  Structure  in  Interconnected  Dynamic  Sys¬ 
tems 

Structure  and  dynamic  behavior  are  two  of  the  most  fundamental  concepts  characterizing 
a  system.  The  interplay  between  these  concepts  has  been  a  central  theme  in  control  as  far 
back  as  Black’s,  Bode’s,  or  Nyquist’s  work  on  feedback  amplifiers  [53,  54,  55],  or  even  as 
early  as  Maxwell’s  analysis  On  Governors  in  1868  [56]. 

The  key  property  driving  such  analyses  is  the  fact  that  an  interconnection  of  systems 
yields  another  system.  This  property  suggests  a  natural  notion  of  structure,  as  the  intercon¬ 
nection  of  systems,  and  focuses  attention  on  understanding  how  interconnections  of  different 
systems  result  in  varieties  of  dynamic  behaviors. 

This  idea  of  structure  as  interconnection  is  not  only  critical  in  the  analysis  of  systems, 
but  it  also  plays  a  key  role  in  system  modeling.  While  black  box  approaches  to  modeling  seek 
to  duplicate  the  input-output  behavior  of  a  system  irrespective  of  structure,  first  principles 
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approaches  to  modeling  use  the  visible  interconnection  structure  of  a  system  to  decompose 
the  system  into  smaller,  fundamental  components.  Constitutive  relations  describing  each 
of  these  components  are  then  applied,  and  interconnection  principles  linking  these  relations 
then  result  in  a  model  structurally  and  dynamically  consistent  with  the  original  system. 

While  such  approaches  have  demonstrated  remarkable  success  in  electrical  and  mechanical 
domains  where  system  components  are  readily  visible  and  physically  separated  from  each 
other  [57,  58,  59],  application  of  such  methods  to  biological,  social,  and  other  domains 
has  been  more  difficult.  One  reason  may  be  that  these  systems  do  not  exhibit  a  natural 
structure  in  the  same  sense  that  previous  applications  have;  while  components  of  electrical 
and  mechanical  systems  are  compartmentalized  and  solid-state,  the  physical  relationship 
among  components  of  these  other  systems  are  often  much  more  fluid  [60,  61].  Perhaps  for 
these  other  domains  different  notions  of  structure  play  the  role  historically  occupied  by  the 
interconnection  of  components. 

This  section  explores  these  ideas  by  characterizing  the  complete  computational  structure 
of  a  system  and  then  contrasting  it  with  three  distinct  partial  structure  representations. 
These  different  representations  include  the  interconnection  of  subsystems  and  the  standard 
idea  of  a  transfer  function  matrix,  but  it  also  includes  a  newer  concept  of  system  structure 
called  signal  structure  that  appears  to  be  especially  useful  for  characterizing  systems  that 
are  difficult  to  compartmentalize.  Precise  relationships  between  these  various  perspectives 
of  system  structure  are  then  provided,  along  with  a  brief  discussion  on  their  implications  for 
various  questions  about  realization  and  approximation. 

4.5.1  Complete  Computational  Structure 

The  complete  computational  structure  of  a  system  characterizes  the  actual  processes  it  uses 
to  sense  properties  of  its  environment,  represent  and  store  variables  internally,  and  affect 
change  externally.  At  the  core  of  these  processes  are  information  retrieval  issues  such  as  the 
encoding,  storage,  and  decoding  of  quantities  that  drive  the  system’s  dynamics.  Different 
mechanisms  for  handling  these  quantities  result  in  different  system  structures. 

Mathematically,  state  equations,  or  their  generalization  as  descriptor  systems  [62,  63,  64] 
are  typically  used  to  describe  these  mechanisms.  Although  there  may  be  many  realizations 
that  describe  the  same  input-output  properties  of  a  particular  system,  its  complete  com¬ 
putational  structure  is  the  architecture  of  the  particular  realization  fundamentally  used  to 
store  state  variables  in  memory  and  transform  system  inputs  to  the  corresponding  outputs. 
In  this  work  we  will  focus  our  attention  on  a  class  of  differential  algebraic  systems  that  are 
equivalent  to  a  set  of  ordinary  differential  equations  in  state  space  form;  we  will  refer  to  such 
equations  as  generalized  state  equations. 

Representing  a  system’s  complete  computational  structure  is  thus  a  question  of  graphi¬ 
cally  representing  the  structure  implied  by  the  equations  that  govern  its  state  evolution.  In 
this  work,  rather  than  focusing  on  the  specific  syntax  of  any  one  particular  graphical  mod¬ 
eling  language,  we  will  draw  from  the  standard  system  theoretic  notions  of  a  block  diagram 
and  a  signal  flow  graph  to  conduct  a  concrete  analysis  between  graphical  representations 
of  a  system  at  various  levels  of  abstraction.  The  complete  computational  structure  of  a 
system,  then,  is  the  description  of  the  system  with  the  most  refined  resolution,  which  we 
will  characterize  as  a  graph  derived  from  a  particular  block  diagram  of  the  generalized  state 
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equations. 

To  make  this  concept  of  structure  precise,  we  begin  by  considering  a  system  G  with 
generalized  state  space  realization 


X  =  f(x,w,u), 

w  =  g(x,w,u),  (44) 

y  =  h(x,w,u). 

Note  that  this  system  is  in  the  form  of  a  differential  algebraic  equation,  although  we  will 

only  consider  systems  with  differentiation  index  zero,  implying  that  (44)  is  always  equivalent 
to  a  standard  ordinary  differential  or  difference  equation  of  the  same  order  [65].  Typically 
we  may  consider  the  system  (44)  to  be  defined  over  continuous  time,  with  1  6  I  >  0,  and 
with  u  G  Mm,  x  G  Rn,  w  G  M1,  y  G  Kp,  and  x  taken  to  mean  dx/dt.  Moreover,  we  restrict 
our  attention  to  those  functions  /,  g  and  h  where  solutions  exist  for  t  >  0.  Nevertheless,  we 
could  also  consider  discrete  time  systems,  with  t  =  0, 1,  2,  3, ...  and  x  in  (44)  taken  to  mean 
x[t  +  1],  or  systems  with  general  input,  state,  auxiliary,  and  output  spaces  U,  X,  W,  or  y , 
respectively.  In  some  situations  these  “spaces”  may  merely  be  sets,  e.g.  X  =  {0, 1}.  In  any 
case,  however,  we  will  take  u  G  Um,  x  G  Xn,  w  G  Wl,  and  y  G  yp  so  that  m,  n,  l  and  p 
characterize  the  dimensions  of  the  input,  state,  auxiliary  and  output  vectors,  respectively. 

Note  that  the  auxiliary  variables,  w,  are  used  to  characterize  intermediate  computation 
in  the  composition  of  functions.  Thus,  for  example,  we  distinguish  between  f(x)  =  x  and 
f(x)  =  2(.5x)  by  computing  the  latter  as  f(w)  =  2w  and  w  =  g(x)  =  .5x.  In  this  way, 
the  auxiliary  variables  serve  to  identify  stages  in  the  computation  of  the  state  space  real¬ 
ization  (44).  Frequently  we  may  not  require  any  auxiliary  variables  in  onr  description  of 
the  system;  indeed  it  is  the  standard  practice  to  eliminate  auxiliary  variables  to  simplify 
the  state  descriptions  of  systems.  Nevertheless,  as  we  discuss  structure,  it  will  be  critical 
to  use  auxiliary  variables  to  distinguish  between  systems  with  dynamically  equivalent,  yet 
structurally  distinct  architectures,  leading  to  the  following  definition. 

Definition  4.  Given  a  system  (44)  >  we  ca M  the  number  of  auxiliary  variables,  l,  the  intricacy 
of  the  realization. 

To  understand  the  structure  of  (44),  we  need  a  notion  of  dependence  of  a  function  on  its 
arguments.  For  example,  the  function  f(x,y,z )  =  xy  —  x  +  z  clearly  depends  on  z,  but  it 
only  depends  on  x  when  y  1  (or  on  y  when  x  0).  Since  “structure”  refers  at  some  level 
to  the  dependence  of  the  system  variables  on  each  other,  it  is  important  that  our  notion  of 
dependence  be  made  clear. 

Definition  5.  A  function  f(w),  from  l-dimensional  domain  W  to  s-dimensional  co-domain 
Z,  is  said  to  depend  on  the  ith  variable,  w if  there  exist  values  of  the  other  s  —  1  variables 
Wj,  j  i,  such  that  f(w)  is  not  constant  over  all  values  of  wy  while  holding  the  others 
variables  fixed.  If  s  =  1,  then  f(w )  depends  on  w  if  it  is  not  constant  over  all  values  ofw. 

Note  that  when  df/dwt  is  well  defined,  the  above  definition  coincides  with  the  partial 
derivative  being  non-zero  for  some  value  of  the  variables  Wj.  Nevertheless,  here  we  allow 
for  non-differentiablc  functions  as  we  explicitly  characterize  one  notion  of  the  structure  of  a 
state  space  realization. 
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Definition  6.  Given  a  system  G  with  realization  (44) >  it s  complete  or  computational  struc¬ 
ture  is  a  weighted  directed  graph  C  with  vertex  set  V{C),  and  edge  set  E(C).  The  vertex  set 
contains  m  +  n  +  l+p  elements,  one  associated  with  the  mechanism  that  produces  each  input, 
state,  auxiliary,  and  output  variable  of  the  system,  and  we  label  the  vertices  accordingly.  In 
particular,  the  vertex  associated  with  the  ith  input  is  labeled  Ui,  1  <  i  <  m,  the  vertex  asso¬ 
ciated  with  the  jth  state  is  labeled  f),  0  <  j  <  n,  the  vertex  associated  with  the  jth  auxiliary 
variable  is  labeled  gv  0  <  j  <  l,  and  the  vertex  associated  with  the  kth  output  is  labeled  h\~, 
1  <  k  <  p.  The  edge  set  contains  an  edge  from  node  i  to  node  j  if  the  function  associated 
with  the  label  of  node  j  depends  on  the  variable  produced  by  node  i.  Moreover,  the  edge  ( i,j ) 
is  then  labeled  (weighted)  with  the  variable  produced  by  node  i. 


So,  for  example,  consider  the  following  continuous  time  system  with  real-valued  variables: 


1  =  fi(xuw,u) 

X2  \  h{x i) 

w  =  g(x2) 

V  =  h(x 2)  . 


(45) 


Its  complete  computational  structure  has  node  set  V(C)  =  {u,  fi,  f2,  g,  h}  and  edge  set 
E(C)  =  {u(u,f1),x1(f1,f1),x1(f1,f2),X2(f2,g),X2(f2,h),w(g,f1)}\  Figure  21  illustrates  the 
graph. 

This  notion  of  the  complete  computational  structure  of  the  system  (44)  corresponds  with 
the  traditional  idea  of  the  structure  of  a  dynamic  system  (see  for  example  [66]),  but  with 
some  important  differences.  First,  this  description  uses  auxiliary  variables  to  keep  track  of 
the  structural  differences  introduced  by  the  composition  of  functions.  This  allows  a  degree 
of  flexibility  in  how  refined  a  view  of  the  computational  structure  one  considers  “complete.” 
Also,  this  description  may  have  some  slight  differences  in  the  labeling  of  nodes  and  edges 
compared  with  various  descriptions  in  the  literature.  These  differences  will  become  important 
as  new  notions  of  partial  structure  are  introduced  in  later  sections.  Here,  we  use  rectangular 
nodes  for  graphs  where  the  nodes  represent  systems ,  and  the  associated  edges  will  represent 
signals.  This  convention  will  bridge  well  between  the  graphical  representation  of  system 
structure  and  the  typical  engineering  block  diagram  of  a  system,  and  it  sometimes  motivates 
the  simplification  in  drawing  edges  as  shown  in  Figure  21b),  since  every  edge  that  leaves  a 
node  represents  the  same  variable  and  carries  the  same  label.  Moreover,  notice  that  nodes 
associated  with  the  mechanisms  that  produce  output  variables  are  terminal,  in  that  no  edge 
ever  leaves  these  nodes,  while  the  nodes  associated  with  input  variables  are  sources,  in  that 
no  edge  ever  arrives  at  these  nodes.  Although  it  is  common  for  engineering  diagrams  to 
explicitly  draw  the  edges  associated  with  output  variables  and  leave  them  “dangling,”  with 
no  explicit  terminal  node,  or  to  eliminate  the  input  nodes  and  simply  depict  the  input  edges- 
also  “dangling,”  our  convention  ensures  that  the  diagram  corresponds  to  a  well  defined  graph, 
with  every  edge  characterized  by  an  ordered  pair  of  nodes.  Note  also  that  state  nodes,  such 
as  f\  in  Figure  21,  may  have  self  loops,  although  auxiliary  nodes  will  not,  and  at  times  it 
will  be  convenient  to  partition  the  vertex  set  into  groups  corresponding  to  the  input,  state, 
auxiliary,  and  output  mechanisms  as  V (C)  =  (VU(C),  VX(C),  VW(C),  Vy(C)j.  Likewise,  we  may 
similarly  partition  the  edge  set  as  necessary. 

We  see,  then,  that  knowing  the  complete  structure  C  of  a  system  is  equivalent  to  knowing 
its  state  space  realization,  along  with  the  composition  structure  with  which  these  functions 
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X, 


(a)  The  complete  computational  structure  C  of  the 
simple  example  specified  by  equation  (45). 


(b)  A  modified  representation  of  the  complete  com¬ 
putational  structure  of  the  system  specified  by 
equation  (45).  Since  the  edges  leaving  a  partic¬ 
ular  node  will  always  represent  the  same  variable, 
they  have  been  combined  to  simplify  the  figure. 


Figure  21:  The  complete  computational  structure  of  the  state  realization  of  a  system  is 
a  graph  demonstrating  the  dependency  among  all  system  variables;  edges  correspond  to 
variables  and  nodes  represent  constitutive  mechanisms  that  produce  each  variable.  System 
outputs  are  understood  to  leave  the  corresponding  terminal  nodes,  and  system  inputs  arrive 
at  the  corresponding  source  nodes. 
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are  represented,  given  by  ( f,g,h ).  We  refer  to  this  structure  as  computational  because  it 
reveals  the  dependencies  among  variables  in  the  particular  representation,  or  basis,  that 
they  are  stored  in  and  retrieved  from  memory.  These  specific,  physical  mechanisms  that 
store  and  retrieve  information,  identified  with  Vx  (C),  are  interconnected  with  devices  that 
transform  variables,  identified  with  VW(C),  and  with  devices  that  interface  with  the  system’s 
external  environment.  These  devices  include  sensors,  identified  with  VU(C),  and  actuators, 
identified  with  Vy(C),  to  implement  the  particular  system  behavior  observed  by  the  outside 
world  through  the  manifest  variables,  u  and  y.  Although  other  technologies  very  well  may 
implement  the  same  observed  behavior  via  a  different  computational  structure  and  a  different 
representation  of  the  hidden  variables,  x  and  w,  C  describes  the  structure  of  the  actual  system 
employing  existing  technologies  as  captured  through  a  particular  state  description.  In  this 
sense,  C  is  the  complete  architecture  of  the  system,  and  often  may  be  interpreted  as  the 
system’s  “physical  layer.”  Importantly,  it  is  often  this  notion  of  structure,  or  a  very  related 
concept,  that  is  meant  when  discussing  the  “structure”  of  a  system,  as  the  next  example 
illustrates. 

Example:  Graph  Dynamical  Systems  As  an  example,  we  examine  the  computational 
structure  of  a  graph  dynamical  system  (GDS).  Graph  dynamical  systems  are  finite  alphabet, 
discrete  time  systems  with  dynamics  defined  in  terms  of  the  structure  of  an  associated 
undirected  graph.  They  have  been  employed  in  the  study  of  various  complex  systems  [67,  68], 
including 

•  Dynamical  process  on  networks: 

—  disease  propagation  over  a  social  contact  graph, 

—  packet  flow  in  cell  phone  communication, 

—  urban  traffic  and  transportation; 

•  Computational  algorithms: 

—  Gauss-Seidel, 

—  gene  annotation  based  on  functional  linkage  networks, 

—  transport  computations  on  irregular  grids; 

•  Computational  paradigms  related  to  distributed  computing. 

Here  we  observe  that  the  computational  structure  of  the  graph  dynamical  system  corresponds 
naturally  with  the  system’s  underlying  graph. 

Given  an  undirected  graph  Q  with  vertex  set  V(Q)  =  {1,2,  ...,n},  a  GDS  associates  with 
each  node  i  a  state  Xi  that  takes  its  values  from  a  specified  finite  set  X .  This  state  is  then 
assigned  a  particular  update  function  that  updates  its  value  according  to  the  values  of  states 
associated  with  nodes  adjacent  to  node  i  on  Q.  Notice  that  this  restriction  on  the  update 
function  suggests  that  the  update  function  for  state  i  depends  on  states  consistent  with  the 
structure  of  Q,  indicating  that  the  system’s  computational  structure  C  should  correspond  to 
the  adjacency  structure  of  Q. 
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(a)  Undirected  graph  Q  defining  the  sequential  (b)  Computational  structure  C  of  the  sequen- 

Graph  Dynamical  System  (46).  tial  Graph  Dynamical  System  (46). 

Figure  22:  Graph  Dynamical  Systems  are  finite  alphabet  discrete  time  systems  with  dynam¬ 
ics  defined  in  terms  of  the  adjacency  structure  of  a  specified  undirected  graph.  Here  we  note 
that  the  computational  structure  of  the  system  reflects  the  structure  of  its  underlying  graph. 


The  distinction  is  made  between  a  GDS  that  updates  all  of  its  states  simultaneously, 
called  a  parallel  GDS ,  and  one  that  updates  its  states  in  a  particular  sequence,  called  a 
sequential  GDS.  The  parallel  GDS  thus  becomes  an  autonomous  dynamical  system,  evolving 
its  state  according  to  its  update  function  from  some  initial  condition.  The  sequential  GDS,  on 
the  other  hand,  can  be  viewed  as  a  controlled  system  that  receives  a  particular  permutation 
of  the  node  set  V(Q)  as  input  and  then  evolves  its  states  accordingly. 

To  illustrate,  consider  the  sequential  GDS  given  by  Q  =  {1,  2,  3, 4}  as  indicated  in  Figure 
22.  Let  U  =  {1, 2,  3, 4}  and  X  =  y  =  {0, 1},  with  update  function  given  by 


where 


Si(x,  n) 


X 1  [t  +  1] 

/i(z[£],u[f]) 

x2[t  +  1] 

f2(x[t},u[t\) 

x3[t  +  1] 

f3(x[t\,u[t\) 

x4[t  +  1]  _ 

f4{x[t],u[t]) 

y[t\ 

=  Xi[t\ 

[  xi 

u  ^  i 

\  (1  +0^(1  +xi_1)(l  TZi+i)  u  =  i 

(46) 


(47) 


and  the  arithmetic  in  (47)  is  taken  modulo  2,  while  that  in  the  subscript  notation  is  taken 
modulo  4  (resulting  in  x$  =  aq,  etc.).  So,  for  example,  the  initial  condition  rc[0]  =  [000  0]T 
with  input  sequence  u[t]  =  1,  2, 3, 4, 1,  2, 3,  4, ...  would  result  in  the  following  periodic  trajec- 
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tory: 


x[l]  =  [1  0  0  0]T 

x[12]  =  [010  0]T 

1 

1 

x[2\  =  [1  0  0  0]T 

x[16]  =  [001  0]T 

1 

1 

x[3]  =  [101  0]T 

x[20]  =  [1  0  0  0]T 

1 

1 

a;  [4]  =  [101  0]T 

x[24]  =  [010  If 

1 

i 

x[8]  =  [0  0  0  1]T 

x[28]  =  [0  0  0  0]T 

The  computational  structure  of  the  system  (46)  follows  immediately  from  the  dependency 
among  variables  characterized  by  equation  (47);  Figure  22  illustrates  C  for  this  system.  No¬ 
tice  that  the  structure  of  Q  is  reflected  in  C,  where  the  undirected  edges  in  Q  have  been 
replaced  by  directed  edges  in  both  directions,  and  self-loops  have  appeared  where  applica¬ 
ble.  Moreover,  notice  that  C  considers  the  explicit  influence  of  the  input  sequence  on  the 
update  computation,  and  explicitly  identifies  the  system  output.  In  this  way,  we  see  that 
the  computational  structure  is  a  reflection  of  the  natural  structure  of  the  sequential  GDS. 

Computational  Structure  of  Linear  Systems  Linear  systems  represent  an  important 
special  case  of  those  described  by  (44).  They  arise  naturally  as  the  linearization  of  sufficiently 
smooth  nonlinear  dynamics  near  an  equilibrium  point  or  limit  cycle,  or  as  the  fundamental 
dynamics  of  systems  engineered  to  behave  linearly  under  nominal  operating  conditions.  In 
either  case,  knowing  the  structure  of  the  relevant  linear  system  is  a  critical  first  step  to 
understanding  that  of  the  underlying  nonlinear  phenomena. 

The  general  state  description  of  a  linear  system  is  given  by 

x  =  Ax  +  Aw  +  Bu , 

w  =  Ax  +  Aw  +  Bu,  (48) 

y  =  Cx  +  Cw  +  Du, 

where  A  G  Knxn,  A  G  Rnxl,  A  €  Rlxn,  A  €  Rlxl,  B  G  Mnxm,  b  G  Rlxm,  C  €  Rpxn,  C  €  Rpxl, 
and  D  G  Mpxm.  Note  that  I  —  A  is  necessarily  invertible,  ensuring  that  the  differentiability 
index  of  the  system  is  zero.  Nevertheless,  the  matrices  are  otherwise  free. 

As  in  the  nonlinear  case,  it  should  be  apparent  that  the  auxiliary  variables  are  superfluous 
in  terms  of  characterizing  the  dynamic  behavior  of  the  system;  this  idea  is  made  precise  in 
the  following  lemma.  Nevertheless,  the  auxiliary  variables  make  a  very  important  difference 
in  terms  of  characterizing  the  system’s  complete  computational  structure,  as  illustrated  by 
the  subsequent  example. 

Lemma  2.  For  any  system  (48)  with  intricacy  l  >  0,  there  exists  a  unique  minimal  intricacy 
realization  (Aa,  B0,C0,  D0)  with  l  =  0  such  that  for  every  solution  (u(t),x(t),w(t),y(t))  of 
(4:8),  (u(t),x(t),y(t))  is  a  solution  of  (A0,  B0,C0,  D0). 

Proof.  The  result  follows  from  the  invertibility  of  (I  —  A).  Solving  for  w  and  substituting 
into  the  equations  of  x  and  y  then  yields  (Aa,  B0,  CQ,  D0).  □ 
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Consider,  for  example,  the  system  (48)  with  state  matrices 
following: 


A 


0  0 
0  0 


0  0  10  10 
0  0  0  1  0  1 


given  by  D 

1  0  ' 

0  1 


0  and  the 


A  = 


ci  0 

0  c2 

0  0 
0  0 
0  0 
0  0 
ai  0 
0  a2 


A 


0  0000000 
0  0000000 
0  Cl  0  0  0  0  0  0 
e2  0000000 
0  0000000 
0  0000000 
0  0000000 
0  0000000 


B 


0  0 
0  0 


bt 


0000  h  000 
00000  b2  00 


c 


0  0 
0  0 


1  0  0  0  0  0  0  0 
0  1  0  0  0  0  0  0 


(49) 


This  system  has  the  complete  computational  structure  C  shown  in  Figure  23a.  Here,  because 
each  auxiliary  variable  is  defined  as  the  simple  product  of  a  coefficient  times  another  variable, 
we  label  the  node  corresponding  to  Wi  in  C  with  the  appropriate  coefficient  rather  than  the 
generic  label,  g*.  Note  that  this  realization  has  an  intricacy  of  /  =  8. 

Suppose,  however,  that  we  eliminate  the  last  six  auxiliary  variables,  leading  to  an  equiv¬ 
alent  realization  with  intricacy  1  =  2.  The  state  matrices  then  become 


A 


a  i  0 

0  a2 


A 


0  ei 

e2  0 


A 


ci  0 

0  c2 


A 


0  0 
0  0 


B 


bi  0 

0  b2 


Bt 


0  0 
0  0 


(50) 


c 


0  0 
0  0 


c 


1  0 
0  1 


with  computational  structure  C  as  shown  in  Figure  23b.  Similarly,  we  can  find  an  equivalent 
realization  with  l  =  0  given  by 


4,= 

Cl\  CiC2 

B0  = 

o 

rO  C 

C0  = 

Ci  0 

e2ci  a2 

1  U  »2  1 

0  c2 

(51) 
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and  all  other  system  matrices  equal  to  zero.  This  realization,  (51),  is  the  minimal  intricacy 
realization  of  both  systems  (49)  and  (50),  and  its  complete  computational  structure  C  is  given 
in  Figure  23c.  The  equivalence  between  these  realizations  is  easily  verified  by  substitution. 

Comparing  the  computational  structures  for  different  realizations  of  the  same  system, 
(49),  (50),  and  (51),  we  note  that  the  intricacy  of  auxiliary  variables  plays  a  critical  role  in 
suppressing  or  revealing  system  structure.  Moreover,  note  that  auxiliary  variables  can  change 
the  nature  of  which  variables  are  manifest  or  hidden.  In  Figure  23  shaded  regions  indicate 
which  variables,  represented  by  edges,  are  hidden;  manifest  variables  leave  shaded  regions 
while  hidden  variables  are  contained  within  them.  Note  that  the  minimal  intricacy  realization 
has  no  internal  manifest  variables,  or,  in  other  words,  it  has  a  single  block  of  hidden  variables 
(Figure  23c).  Meanwhile,  both  w\  and  W2  are  manifest  in  the  other  realizations  (Figures  23a, 
23b)  since  w i  =  y\  and  W2  =  1/2,  indicated  by  the  “1”  at  their  respective  terminal  nodes. 
This  yields  two  distinct  blocks  of  hidden  variables,  in  either  case,  revealing  the  role  intricacy 
of  a  realization  can  play  characterizing  its  structure. 

The  complete  computational  structure  of  a  system  is  thus  a  graphical  representation  of 
the  dependency  among  input,  state,  auxiliary,  and  output  variables  that  is  in  direct,  one-to- 
one  correspondence  with  the  system’s  state  realization,  generalized  to  explicitly  account  for 
composition  intricacy.  All  structural  and  behavioral  information  is  fully  represented  by  this 
description  of  a  system.  Nevertheless,  this  representation  of  the  system  can  also  be  unwieldy 
for  large  systems  with  intricate  structure. 

4.5.2  Partial  Structure  Representations 

Complex  systems  are  often  characterized  by  intricate  computational  structure  and  compli¬ 
cated  dynamic  behavior.  State  descriptions  and  their  corresponding  complete  computational 
structures  accurately  capture  both  the  system’s  structural  and  dynamic  complexity,  never¬ 
theless  these  descriptions  themselves  can  be  too  complicated  to  convey  an  efficient  under¬ 
standing  of  the  nature  of  the  system.  Simplified  representations  are  then  desirable. 

One  way  to  simplify  the  representation  of  a  system  is  to  restrict  the  structural  information 
of  the  representation  while  maintaining  a  complete  description  of  the  system’s  dynamics.  The 
most  extreme  example  of  this  type  of  simplified  representation  is  the  transfer  function  of  a 
single-input  single-output  linear  time  invariant  (LTI)  system.  A  transfer  function  completely 
specifies  the  system’s  input-output  dynamics  without  retaining  any  information  about  the 
computational  structure  of  the  system.  For  example,  consider  the  nth  order  LTI  single-input 
single-output  system  given  by  (A,  b ,  c,  d).  It  is  well  known  that  although  the  state  description 
of  the  system  completely  specifies  the  transfer  function,  G(s)  =  c(sl  —  A)~lb  +  d,  the 
transfer  function  G(s)  has  an  infinite  variety  of  state  realizations,  and  hence  computational 
structures,  that  all  characterize  the  same  input-output  behavior.  That  is,  the  structural 
information  in  any  state  realization  of  the  system  is  completely  removed  in  the  transfer 
function  representation  of  the  system,  even  though  the  dynamic  (or  behavioral)  information 
about  the  system  is  preserved. 

We  use  this  power  of  a  transfer  function  to  obfuscate  structural  information  to  develop 
three  distinct  partial-structure  representations  of  an  LTI  system:  subsystem  structure,  signal 
structure,  and  the  zero  pattern  of  a  (multiple  input,  multiple  output)  system’s  transfer 
function  matrix.  Later  we  will  show  how  each  of  these  representations  has  more  structural 
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I 


(a)  The  complete  computational  structure  C  of  the 
linear  system  given  by  the  state  matrices  (49)  with 
intricacy  l  =  8. 


l_J 

1*71 

/ 


(b)  The  complete  computational  structure  C  of  the 
equivalent  linear  system  with  intricacy  1  =  2,  spec¬ 
ified  by  the  realization  (50). 


(c)  The  complete  computational  structure  C  of  the 
minimal  intricacy  (Z  =  0)  realization  for  both  sys¬ 
tems  above,  characterized  by  (51). 

Figure  23:  The  complete  computational  structure  of  a  linear  system  characterized  by  realiza¬ 
tions  of  differing  intricacies.  Edges  within  shaded  regions  represent  hidden  variables,  while 
those  outside  shaded  regions  are  manifest  variables. 
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information  than  its  successor,  and  we  will  precisely  characterize  the  relationships  among 
them. 

Subsystem  Structure  One  of  the  most  natural  ways  to  reduce  the  structural  information 
in  a  system’s  representation  is  to  partition  the  nodes  of  its  computational  structure  into  sub¬ 
systems,  then  replace  these  subsystems  with  their  associated  transfer  function.  Each  transfer 
function  obfuscates  the  structure  of  its  associated  subsystem,  and  the  remaining  (partial) 
structural  information  in  the  system  is  the  interconnection  between  transfer  functions. 

Subsystem  structure  refers  to  the  appropriate  decomposition  of  a  system  into  constituent 
subsystems  and  the  interconnection  structure  between  these  subsystems.  Abstractly,  it  is  the 
condensation  graph  of  the  complete  computational  structure  graph,  C,  taken  with  respect  to 
a  particular  partition  of  C  that  identifies  subsystems  in  the  system.  Such  abstractions  have 
been  used  in  various  ways  to  simplify  the  structural  descriptions  of  complex  systems  [66,  69], 
for  example  by  “condensing”  strongly  connected  components  or  other  groups  of  vertices  of 
a  graph  into  single  nodes.  Nevertheless,  in  this  work  we  define  a  particular  condensation 
graph  as  the  subsystem  structure  of  the  system.  We  begin  by  characterizing  the  partitions 
of  C  that  identify  subsystems. 

Definition  7.  Given  a  system  G  with  realization  (48)  and  associated  computational  structure 
C,  we  say  a  partition  ofV(C)  is  admissible  if  every  edge  in  E(C )  between  components  of  the 
partition  represents  a  variable  that  is  manifest,  not  hidden. 

For  example,  considering  the  system  (51)  with  V ( C )  =  {wi,  /i,  c i,  c2,  f2,  u2}.  We  see  that 
the  partition  {(«i),  (/i,  Ci,  c2, /2),  (u2)}  is  admissible  since  the  only  edges  between  compo¬ 
nents  are  Ui(mi,/i)  and  u2(u2,f2),  representing  the  manifest  variables  and  u2.  Notice 
that  the  shading  in  Figure  23c  is  consistent  with  this  admissible  partition.  Alternatively, 
the  partition  {(«i),  (/i,  Ci),  (c2,  /2),  (w2)}  is  not  admissible  for  (51),  since  the  edges  /2) 

and  x2(/2,  fi )  extend  between  components  of  the  partition  but  represent  variables  X\  and  x2 
that  are  hidden,  not  manifest. 

Although  sometimes  any  aggregation,  or  set  of  fundamental  computational  mechanisms 
represented  by  vertices  in  C ,  may  be  considered  a  valid  subsystem,  in  this  work  a  subsystem 
has  a  specific  meaning.  In  particular,  the  variables  that  interconnect  subsystems  must  be 
manifest,  and  thus  subsystems  are  identified  by  the  components  of  admissible  partitions 
of  V[C).  We  adopt  this  convention  to  1)  enable  the  distinction  between  real  subsystems 
vs.  merely  arbitrary  aggregations  of  the  components  of  a  system,  and  2)  ensure  that  the 
actual  subsystem  architecture  of  a  particular  system  is  adequately  reflected  in  the  system’s 
computational  structure  and  associated  realization,  thereby  ensuring  that  such  realization  is 
complete. 

Definition  8.  Given  a  system  G  with  realization  (48)  and  associated  computational  structure 
C,  the  system’s  subsystem  structure  is  a  condensation  graph  S  of  C  with  vertex  set  V(S)  and 
edge  set  E(S)  given  by: 

•  F(iS)  =  {Ai,...^}  are  the  elements  of  an  admissible  partition  of  V (C)  of  maximal 
cardinality,  and 
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(a)  The  subsystem  structure  S  of  the  linear  sys¬ 
tem  given  by  the  state  matrices  (49)  with  intri¬ 
cacy  1=8.  Note  that  the  subsystem  structure  of 
the  equivalent,  but  less  intricate  system  given  by 
equation  (50)  is  exactly  the  same,  indicated  by  the 
shaded  regions  in  Figures  23a  and  23b. 


(b)  The  subsystem  partial  structure  S  of  the  dy¬ 
namically  equivalent,  minimally  intricate  linear 
system,  specified  by  the  realization  (51)  corre¬ 
sponding  to  Figure  23c. 


Figure  24:  The  subsystem  structure  of  a  linear  system  partitions  vertices  of  the  complete 
computational  structure  and  condenses  admissible  groups  of  nodes  into  single  subsystem 
nodes,  shown  in  brown;  nodes  shown  in  green  are  those  that  correspond  directly  to  unaggre¬ 
gated  vertices  of  the  complete  computational  structure,  which  will  always  be  associated  with 
mechanisms  that  generate  manifest  variables.  Here,  the  subsystem  structure  corresponds  to 
the  interconnection  of  shaded  regions  in  Figure  23.  Comparing  Figures  24a  and  24b,  we  note 
that  intricacy  in  a  realization  may  be  necessary  to  characterize  meaningful  subsystems  and 
yield  nontrivial  subsystem  structure. 


•  E(S)  has  an  edge  ( Si,Sj )  if  E{C)  has  an  edge  from  some  component  of  Si  to  some 
component  of  Sj . 

We  label  the  nodes  of  V (S)  with  the  transfer  function  of  the  associated  subsystem,  which  we 
also  denote  Sit  and  the  edges  of  E(S)  with  the  associated  variable  from  E(C). 

Note  that,  like  C,  the  subsystem  structure  S  is  a  graph  with  vertices  that  represent 
systems  and  edges  that  represent  signals ,  or  system  variables.  For  example,  Figure  24a 
illustrates  the  subsystem  structure  for  both  systems  (49)  and  (50),  showm  in  Figures  23a  and 
23b.  Note  that  the  subsystem  structure  of  these  systems’  minimally  intricate  realization, 
(51),  is  quite  different,  with  a  single  block  rather  than  two  blocks  interconnected  in  feedback, 
as  shown  in  Figure  24b.  This  illustrates  the  necessity  of  auxiliary  variables  to  adequately 
describe  the  complete  system  structure. 

Lemma  3.  The  subsystem  structure  S  of  a  system  G,  with  complete  computational  structure 
C,  is  unique. 
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Proof.  We  prove  by  contradiction.  Suppose  the  subsystem  structure  S  of  G  is  not  unique. 
Then  there  are  at  least  two  distinct  subsystem  structures  of  G,  which  we  will  label  51  and 
S2.  This  implies  there  are  two  admissible  partitions  of  V(C),  given  by  ^(iS1)  and  I/(iS2), 
such  that  ^(iS1)  ^  V"(52)  and  with  equal  cardinality,  q.  Note  that  by  definition,  q  must  be 
the  maximal  cardinality  of  any  admisiblc  partition  of  V(C).  To  obtain  a  contradiction,  we 
will  construct  another  admissible  partition,  14(<S3),  such  that  |Id(<S3)|  >  q. 

Consider  the  following  partition  of  V(C)  that  is  a  refinement  of  both  C(51)  and  C(52): 

C(53)  =  {S3\S3  72(/}]S3  =  Si  n  Sj,  Si  e  Vis1),  Sj  e  C(52)}. 

Since  y(5x)  ^  C(52),  then \V(S3)\  >  q,  since  the  cardinality  of  the  refinement  must  then  be 
greater  than  that  of  C(51)  or  V(S2).  Moreover,  note  that  the  partition  C(53)  is  admissible, 
since  every  edge  of  C  between  vertices  associated  with  distinct  components  of  Y(<S3)  corre¬ 
sponds  with  an  edge  of  either  <S'  or  S2,  which  are  admissible.  Thus,  C(53)  is  an  admissible 
partition  of  V ( C )  with  cardinality  greater  than  q,  which  contradicts  the  assumption  that  <S' 
and  S2  are  both  subsystem  structures  of  G.  □ 

The  subsystem  structure  of  a  system  reveals  the  way  natural  subsystems  are  intercon¬ 
nected,  and  it  can  be  represented  in  other  ways  besides  (but  equivalent  to)  specifying  S. 
For  example,  one  common  way  to  identify  this  kind  of  subsystem  architecture  is  to  write 
the  system  as  the  linear  fractional  transformation  (LFT)  with  a  block  diagonal  “subsystem” 
component  and  a  static  “interconnection”  component  (see  [70]  for  background  on  the  LFT). 
For  example,  the  system  in  Figure  24a  can  be  equivalently  represented  by  the  feedback  in¬ 
terconnection  of  a  static  system  !V:UxW— >-Yx(Ux  W)  and  a  block-diagonal  dynamic 
system  S:UxWgW  given  by 

0  0  10' 

0  0  0  1 
10  0  0 
0  0  0  1  ’ 

0  0  10 
0  1  0  0  _ 

where 

S\  :  Ul  — *  ui\  S2  :  Wl  — >  w2  (53) 

W2  \  L  U2 

In  general,  the  LFT  associated  with  S  will  have  the  form 


where  q  is  the  number  of  distinct  subsystems,  and  L  and  K  are  each  binary  matrices  of  the 
appropriate  dimension  (see  Figure  25).  Note  that  if  additional  output  variables  are  present, 
besides  the  manifest  variables  used  to  interconnect  subsystems,  then  the  structure  of  N 
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Figure  25:  The  subsystem  structure  of  a  system  can  always  be  represented  by  the  lower 
linear  fractional  transformation  of  the  static  interconnection  matrix  N  with  a  block  diagonal 
transfer  function  matrix  S.  Note  that  7r(u,w)  represents  a  permutation  of  a  subset  of  the 
variables  in  the  vector  inputs,  u,  and  manifest  auxiliary  variables,  w,  possibly  with  repetition 
of  some  variables  if  necessary. 


and  S  above  extend  naturally.  In  any  event,  however,  N  is  static  and  L  and  K  are  binary 
matrices. 

The  subsystem  structure  of  any  system  is  well  defined  although  it  may  be  trivial  (a 
single  internal  block)  if  the  system  does  not  decompose  naturally  into  an  interconnection 
of  subsystems,  such  as  (51)  in  Figure  23c  and  Figure  24b.  Note  that  S  always  identifies 
the  most  refined  subsystem  structure  possible,  and  for  systems  with  many  interconnected 
subsystems,  coarser  representations  may  be  obtained  by  aggregating  various  subsystems  to¬ 
gether  and  constructing  the  resulting  condensation  graph.  These  coarser  representations 
effectively  absorb  some  interconnection  variables  and  their  associated  edges  into  the  aggre¬ 
gated  components,  suggesting  that  such  representations  are  the  subsystem  structure  for  less 
intricate  realizations  of  the  system,  where  some  of  the  manifest  auxiliary  variables  are  re¬ 
moved,  or  at  least  left  hidden.  The  subsystem  structure  is  thus  a  natural  partial  description 
of  system  structure  when  the  system  can  be  decomposed  into  the  interconnection  of  specific 
subsystems. 

Signal  Structure  Another  very  natural  way  to  partially  describe  the  structure  of  a  system 
is  to  characterize  the  direct  causal  dependence  among  each  of  its  manifest  variables;  we  will 
refer  to  this  notion  as  the  signal  structure.  This  description  of  the  structure  of  a  system 
makes  no  attempt  to  cluster,  or  partition,  the  actual  internal  system  states.  As  a  result,  it 
offers  no  information  about  the  internal  interconnection  of  subsystems,  and  signal  structure 
can  therefore  be  a  very  different  description  of  system  structure  than  subsystem  structure. 

Given  a  generalized  linear  system  (48)  with  complete  computational  structure  C,  we  char¬ 
acterize  its  signal  structure  by  considering  its  minimal  intricacy  realization  (A0,  Ba,  C0,  D0 ). 
We  assume  without  loss  of  generality  that  the  outputs  y  are  ordered  in  such  a  way  that  Co 


APPROVED  FOR  PUBLIC  RELEASE:  DISTRIBUTION  UNLIMITED 

52 


can  be  partitioned 


Cn 


C\  1  C 12 
Cr2l  C22 


where  Cn  G  MpiXpi  is  invertible,  with  p\  equal  to  the  rank  of  CQ\  C\2  G  ]g>pix(n-pi).  (j,n  g 
j^(p-pi)xpi.  anc[  £<22  e  ]g(p-pi)x(n-pi)_  Note  that  if  the  outputs  of  the  minimal  intricacy 
realization  do  not  result  in  C\  \  being  invertible,  then  it  is  possible  to  reorder  them  so  the 
first  pi  outputs  correspond  to  independent  rows  of  C0;  the  states  can  then  be  renumbered 
so  that  C\  1  is  invertible.  One  can  show  that  such  reordering  of  the  outputs  and  states  of 
the  minimal  intricacy  realization  only  affects  the  ordering  of  the  states  and  outputs  of  the 
original  system;  the  graphical  relationship  of  the  computational  structure  is  preserved. 

The  direct  causal  dependence  among  manifest  variables  is  then  revealed  as  follows.  First, 
consider  the  state  transformation  z  =  Tx  given  by 


Cn  C12 

0  / 


This  transformation  yields  a  system  of  the  form 


Zi 

.  ^2  . 

- 

yi 

- 

.  2/2 . 

Mi 

Mi 

I 

C2 


M2 

M2 


Z] 

+ 

'  Bi  ' 

.  Z2  . 

.  B  -2  _ 

u 


0 

0 


Z\ 

+ 

'  D1  ' 

.  Z2  . 

d2 

u 


(55) 


(56) 


where  z  G  M™,  u  G  Mm,  y\  G  MP1,  and  y2  G  W~P1.  To  simplify  the  exposition  we  will  abuse 
notation  and  refer  to  the  above  system  as  (A,  B,C,  D),  since  there  is  little  opportunity  to 
confuse  these  matrices  with  those  of  the  original  system  given  in  (48).  In  fact,  the  system 
(56)  is  simply  a  change  of  coordinates  of  the  minimal  intricacy  realization  (A0,  B0,  CQ,  D0), 
possibly  with  a  reordering  of  the  output  and  state  variables.  The  direct  causal  depen¬ 
dence  among  manifest  variables  is  then  revealed  by  the  dynamical  structure  function  of 

(MB, [I  0  ],D1). 

The  dynamical  structure  function  of  a  class  of  linear  systems  was  defined  in  [71]  and 
discussed  in  [72,  73,  74,  75,  76,  77].  This  representation  of  a  linear  system  describes  the  direct 
causal  dependence  among  a  subset  of  state  variables,  and  it  will  extend  to  characterize  signal 
structure  for  the  system  in  (56).  We  repeat  and  extend  the  derivation  here  to  demonstrate 
its  applicability  to  the  system  (56).  Taking  Laplace  transforms  and  assuming  zero  initial 
conditions  yields  the  following  relationships 


sZ\ 

'  An 

A\2 

'  Z 1  " 

'  Bi 

sZ2 

~ 

A2\ 

a22 

Z2 

+ 

b2 

(57) 


where  Z(s)  denotes  the  Laplace  transform  of  z(t),  etc.  Solving  for  Z2  in  the  second  equation 
and  substituting  into  the  first  then  yields 


sZ  1  =  W(s)Zi  +  V(s)U 


(58) 
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where  W(s)  =  [An  +  Ai2(sl  -  A22 )  1A21\  and  V(s)  =  [B1  +  A12(sl  -  A22)  lB2}.  Let  D(s) 
be  the  matrix  of  the  diagonal  entries  of  W(s),  yielding 

Zx  =  Q(s)Z1  +  P(s)U  (59) 

where  Q(s)  =  (si  —  Z)(s))_1(hL(s)  —  D(s))  and  P(s)  =  (si  —  D(s))~1V(s).  From  (56)  we 
note  that  Z\  =  Y\  —  DiU,  which,  substituting  into  (59),  then  yields: 

'  Y1 

. 

The  matrices  [  Q(s)T  Cf  ]T  and  [  (P(s)  +  (I  —  Q(s))Di)T  Dj  ]T  we  refer  to  as  Q 
and  P,  respectively.  The  matrices  (Q(s),  P(s))  are  called  the  dynamical  structure  function 
of  the  system  (56),  and  they  characterize  the  dependency  graph  among  manifest  variables 
as  indicated  in  Equation  (60).  We  note  a  few  characteristics  of  (Q(s),  P(s))  that  give  them 
the  interpretation  of  system  structure,  namely: 

•  Q(s)  is  a  square  matrix  of  strictly  proper  real  rational  functions  of  the  Laplace  variable, 
s,  with  zeros  on  the  diagonal.  Thus,  if  each  entry  of  y\  were  the  node  of  a  graph,  Qij(s) 
would  represent  the  weight  of  a  directed  edge  from  node  j  to  node  i;  the  fact  Qij(s)  is 
proper  preserves  the  meaning  of  the  directed  edge  as  a  causal  dependency  of  y^  on  yr 

•  Similarly,  the  entries  of  the  matrix  [P(s)  +  (I  —  Q(s))-Di]  carry  the  interpretation  of 
causal  weights  characterizing  the  dependency  of  entires  of  y\  on  the  m  inputs,  u.  Note 
that  when  D\  =  0,  this  matrix  reduces  to  P(s),  which  has  strictly  proper  entries. 

This  leads  naturally  to  the  definition  of  signal  structure. 

Definition  9.  The  signal  structure  of  a  system  G,  with  realization  (48)  and  equivalent 
realization  (56),  and  with  dynamical  structure  function  (Q(s),  P(s))  characterized  by  (59), 
is  a  graph  W,  with  a  vertex  set  V(W)  and  edge  set  E(W)  given  by: 

•  V(W)  =  {ui ,  ...,um,yn, yiPl,  y2\, . . . ,  y2p2 } ,  each  representing  a  manifest  signal  of  the 
system,  and 

•  E(W)  has  an  edge  from  Ui  to  y\v  ut  to  y2j,  y u  to  yij  or  yu  to  y2j  if  the  associated  entry 
in  [P(s)  +  (I  —  Q(s))Di],  D2,  Q(s),  or  C2i  (as  given  in  Equations  (59)  and  (60))  is 
nonzero,  respectively. 

We  label  the  nodes  ofV(W )  with  the  name  of  the  associated  variable,  and  the  edges  of  E(W) 
with  the  associated  transfer  function  entry  from  Equation  (60). 

Signal  structure  is  fundamentally  a  different  type  of  graph  than  either  the  computational 
or  subsystem  structure  of  a  system  because,  unlike  these  other  graphs,  vertices  of  a  system’s 
signal  structure  represent  signals  rather  than  systems.  Likewise,  the  edges  of  W  represent 
systems  instead  of  signals,  as  opposed  to  C  or  5.  We  highlight  these  differences  by  using 
circular  nodes  in  W,  in  contrast  to  the  square  nodes  in  C  or  S.  The  next  example  illustrates  a 
system  without  notable  subsystem  structure  and  no  apparent  structural  motif  in  its  complete 
computational  structure;  nevertheless,  it  reveals  a  simple  and  elegant  ring  structure  in  the 
weak  sense. 


Q(s) 

C2 


w  + 


P(s)  +  (I-Q(s))D1 
D2 


u 


(60) 
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Figure  26:  The  complete  computational  structure  C,  of  the  system  (A0,  Ba,  CQ,  Da)  given  by 
(61).  Here,  edge  labels,  u  and  x,  and  self  loops  on  each  node  fi  have  been  omitted  to  avoid 
the  resulting  visual  complexity.  Edges  associated  with  variables  x,  which  are  not  manifest, 
are  entirely  contained  within  the  shaded  region  (which  also  corresponds  to  the  strong  partial 
condensation  shown  in  Figure  27a). 


Example  3.  Ring  Signal  Structure.  Systems  with  no  apparent  structure  in  any  other  sense 
can  nevertheless  be  very  structured  in  the  weak  sense.  Consider  the  minimally  intricate  linear 
system,  specified  by  the  state-space  realization  (Aa,  B0,C0,  D0),  where 
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(61) 


We  compute  the  signal  structure  by  employing  a  change  of  coordinates  on  the  state  variables 
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(a)  The  subsystem  structure,  S,  of  the  system 
shown  in  Figures  26  and  27b;  the  system  has  no 
interconnection  structure  between  subsystems  be¬ 
cause  the  system  is  composed  of  only  a  single  sub¬ 
system;  it  does  not  meaningfully  decompose  into 
smaller  subsystems. 


(b)  The  signal  structure,  W,  of  the  system  shown 
in  Figures  26  and  27a.  Note  that,  in  contrast  with 
C  and  S,  vertices  of  W  represent  manifest  signals 
(characterized  by  round  nodes),  while  edges  repre¬ 
sent  systems. 

Figure  27:  Distinct  notions  of  partial  structure  for  the  system  (Aa,  B0,C0,  Da)  given  by 
Equation  (61).  Observe  that  while  the  system  is  unstructured  in  the  strong  sense,  it  is 
actually  quite  structured  in  the  weak  sense.  Moreover,  although  the  subsystem  structure  is 
visible  in  the  complete  computational  structure  (Figure  26)  as  a  natural  condensation  (note 
the  shaded  region),  the  signal  structure  is  not  readily  apparent. 
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to  find,  an  equivalent  realization  of  the  form  (56).  The  transformation 


T  = 


-1  4 

-12  21 
0  2 
0  0 
0  0 
0  0 


results  in  the  following  realization 


A  =  TA0T~l 


-2 

0 

0 

1 

0 

0 


B 


TBn  = 


2  0  0 
0  3  0 
0  0  6 
0  0  0 
0  0  0 
0  0  0 


-1  -2  -1  1 
0  -11  -5  6 

0  -2-10 
0  10  0 

0  0  10 

0  0  0  1 


0  0  0  0  -3 

-3  0-10  0 

0-40  5  0 

0  0-40  0 

2  0  0  -5  0 

0  10  0-1 


c  =  CoT-1  =[I3  0  ]  , 


D  =  D0=  [0], 


The  dynamical  structure  function  of  the  system,  ( Q,P ),  then  becomes 


Q 


o 


s2+7s+12 

0 


0 

0 

10 

•s2+9s+20 


s2+3s+2 

o 

0 


(62) 


(63) 


P  = 


2 

s+2 


0 

0 


0 

3 

s+3 

0 


0 

0 

6 

s+4  . 


(64) 


which  yields  the  signal  structure,  W,  as  shown  in  Figure  27b.  Notice  that  although  the  com¬ 
plete  computational  structure  and  subsystem  structure  do  not  characterize  any  meaningful 
interconnection  patterns,  the  system  is  nevertheless  structured  and  organized  in  a  very  con¬ 
crete  sense.  In  particular,  the  outputs  ij\ ,  y2,  and  y3  form  a  cyclic  dependency,  and  each 
causally  depends  on  a  single  input  Ui,i  =  1,2,3,,  respectively. 
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Zero  Pattern  of  the  Transfer  Function  Matrix  The  weakest  notion  of  structure  ex¬ 
hibited  by  a  system  is  the  pattern  of  zeros  portrayed  in  its  transfer  function  matrix,  where 
“zero”  refers  to  the  value  of  the  particular  transfer  function  element,  not  a  transmission 
zero  of  the  system.  Like  signal  structure,  this  type  of  structure  is  particularly  meaning¬ 
ful  for  multiple-input  multiple-output  systems,  and,  like  signal  structure,  the  corresponding 
graphical  representation  reflects  the  dependance  of  system  output  variables  on  system  input 
variables.  Thus,  nodes  of  the  graph  will  be  signals,  represented  by  circular  nodes,  and  the 
edges  of  the  graph  will  represent  systems,  labeled  with  the  corresponding  transfer  function 
element;  a  zero  element  thus  corresponds  to  the  absence  of  an  edge  between  the  associated 
system  input  and  output.  Formally,  we  have  the  following  definition. 

Definition  10.  The  zero  pattern  of  the  transfer  function  of  a  system  G  is  a  graph  Z,  with 
a  vertex  set  V(Z)  and  edge  set  E(Z)  given  by: 

•  V(Z)  =  {«!,  ...,yp},  each  representing  a  manifest  signal  of  the  system,  and 

•  E(Z)  has  an  edge  from  Uj  to  yj  if  Gji  is  nonzero. 

We  label  the  nodes  ofV(Z)  with  the  name  of  the  associated  variable,  and  the  edges  of  E(Z) 
with  the  associated  element  from  the  transfer  function  G(s). 

LInlike  signal  structure,  note  that  the  zero  pattern  of  the  transfer  function  matrix  de¬ 
scribes  the  closed-loop  dependency  of  an  output  variable  on  a  particular  input  variable,  not 
its  direct  dependence.  As  a  result,  the  graph  is  necessarily  bipartite,  and  all  edges  will  begin 
at  an  input  node  and  terminate  at  an  output  node;  no  edges  will  illustrate  dependencies 
between  output  variables.  For  example,  the  zero  pattern  of  the  transfer  function  for  the 
system  in  Example  3,  is  shown  in  Figure  28,  with  transfer  function  G(s)  = 

ni  (s)  — 90(s  +  4)  -18(s3  +  12s2+47s  +  60) 

d(s)  d(s)  d(s) 

-2(s3  +  10s2+29s  +  20)  3(s5+16s4  +  97s3+274s2+352s  +  160)  18(s  +  5) 

- dJT) -  - d{T) -  d(s)  ’  ''DO; 

-20(s  +  l)  30(s3+7s2  +  14s+8)  ra2(s) 

d(s)  d(s)  d(s) 

where  d(s)  =  s6  +  19s5  +  145s4  +  565s3  +  1170s2  +  1220s  +  450,  m(s)  =  2(s5  +  17s4  +  Ills3  + 
343s2  +  488s +  240),  and  n2(s)  =  6(s5  +  15s4  +  85s3  +  225s2  +  274s +  120).  Although  the  direct 
dependence,  given  by  the  signal  structure,  is  cyclic  (Figure  27b),  there  is  a  path  from  every 
input  to  every  output  that  does  not  cancel.  Thus,  the  zero  pattern  of  the  transfer  function 
is  fully  connected,  corresponding  to  the  full  transfer  function  matrix  in  Equation  (65). 

It  is  important  to  understand  that  the  zero  pattern  of  a  transfer  function  does  not  neces¬ 
sarily  describe  the  flow  of  information  between  inputs  and  outputs.  The  presence  of  a  zero 
in  the  (i,  j)th  location  simply  indicates  that  the  input-output  response  of  the  system  results 
in  the  ith  output  having  no  dependence  on  the  jth  input.  Such  an  effect,  however,  could  be 
the  result  of  certain  internal  cancellations  and  does  not  suggest,  for  example,  that  there  is 
no  path  in  the  complete  computational  structure  from  the  jth  input  to  the  ith  output.  Thus, 
for  example,  a  diagonal  transfer  function  matrix  does  not  imply  the  system  is  decoupled. 

The  next  example  demonstrates  that  a  fully  coupled  system  may,  nevertheless,  have  a  di¬ 
agonal  transfer  function,  even  when  the  system  is  minimal  in  both  intricacy  and  order.  That 
is,  the  internal  cancellations  necessary  to  generate  the  diagonal  zero  pattern  in  the  transfer 
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Figure  28:  The  zero  pattern  of  the  transfer  function  for  the  system  in  Example  3.  Note  that 
vertices  of  the  zero  pattern  are  manifest  signals,  distinguished  in  this  work  by  circular  nodes 
similar  to  those  of  signal  structure.  Nevertheless,  this  system  has  a  fully  connected  (and 
thus  unstructured)  zero  pattern,  while  its  signal  structure  (shown  in  Figure  27b)  exhibits  a 
definite  ring  structure. 


function  while  being  fully  coupled  do  not  result  in  a  loss  of  controllability  or  observability. 
Thus,  the  zero  pattern  of  a  transfer  function  only  describes  the  input-output  dependencies 
of  the  system  and  does  not  imply  anything  about  the  internal  information  flow  or  commu¬ 
nication  structure.  This  fact  is  especially  important  in  the  context  of  decentralized  control, 
where  the  focus  is  often  on  shaping  the  zero  pattern  of  a  controller’s  transfer  function  given 
a  particular  zero  pattern  in  the  transfer  function  of  the  system  to  be  controlled. 


Example  4.  Diagonal  Transfer  Function  ^  Decoupled.  Consider  a  system,  G,  with  the 
following  minimal  intricacy  realization,  (Aa,  B0,C0,  D0): 


x\ 

'  -5  1 

x-i 

to 

Ui 

X2 

2  -4 

X2 

+ 

1 

1 

h-L 

U2 

yi 

1  1 ' 

X\ 

V2 

1 

1 

to 

X2 

(66) 


It  is  easy  to  see  from  (A0,  B0 ,  CQ,  D0)  that  this  system  has  a  fully  connected  complete  com¬ 
putational  structure.  Moreover,  one  can  easily  check  that  the  realization  is  controllable  and 
observable,  and  thus  of  minimal  order.  Nevertheless,  its  transfer  function  is  diagonal,  given 
by 


G(s ) 


6 

s+3 


o 


0 

-6 

s+6  . 


(67) 


4.5.3  Relationships  Among  Representations  of  Structure 

In  this  section,  we  explore  the  relationships  between  the  four  representations  of  structure  de¬ 
fined  above.  What  we  will  find  is  that  some  representations  of  structure  are  more  informative 
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than  others.  Next,  we  will  discus  a  specific  class  of  systems  in  which  the  four  representations 
of  structure  can  be  ranked  by  information  content  (with  one  representation  encapsulating 
all  the  information  contained  in  another  representation  of  structure).  Outside  this  class  of 
systems,  however,  we  will  demonstrate  that  signal  structure  and  subsystem  structure  have 
fundamental  differences  in  addition  to  those  arising  from  the  graphical  conventions  of  cir¬ 
cular  versus  square  nodes,  etc.  Signal  and  subsystem  structure  provide  alternative  ways  of 
discussing  a  system’s  structure  without  requiring  the  full  detail  of  a  state-space  realization 
or  the  abstraction  imposed  by  a  transfer  function.  Understanding  the  relationships  between 
these  representations  then  enables  the  study  of  new  kinds  of  research  problems  that  deal 
with  the  realization,  reconstruction,  and  approximation  of  system  structure  (where  struc¬ 
ture  can  refer  to  the  four  representations  defined  so  far  or  other  representations  of  structure 
not  mentioned  in  this  work).  The  final  section  gives  a  discussion  of  such  future  research. 

Different  representations  of  structure  contain  different  kinds  of  structural  information 
about  the  system.  For  example,  the  complete  computational  structure  details  the  struc¬ 
tural  dependencies  among  fundamental  units  of  computation.  Using  complete  computational 
structure  to  model  system  structure  requires  knowledge  of  the  parameters  associated  with 
each  fundamental  unit  of  computation.  Partial  representations  of  structure  such  as  signal 
and  subsystem  structure  do  not  require  knowledge  of  such  details  in  their  description.  Specif¬ 
ically,  subsystem  structure  essentially  aggregates  units  of  computation  to  form  subsystems 
and  models  the  closed-loop  transfer  function  of  each  subsystem.  Signal  structure  models  the 
SISO  transfer  functions  describing  direct  causal  dependencies  between  outputs  and  inputs 
of  some  of  the  fundamental  units  of  computation  that  happen  to  be  manifest.  Zero  pattern 
structure  models  the  closed-loop  dependencies  of  system  outputs  on  inputs.  Thus,  complete 
computational  structure  appears  to  be  the  most  demanding  or  information-rich  description 
of  system  structure.  Indeed,  this  intuition  is  made  precise  with  the  following  result: 


Theorem  1.  Suppose  a  complete  computational  structure  has  minimal  intricacy  realization 
( A0,B0,C0,D0 )  with 

r,  _  Cu  C\2 

U21  022_ 

and  C\  1  invertible.  Then  the  complete  computational  structure  specifies  a  unique  subsystem, 
signal,  and  zero  pattern  structure. 


Proof.  Let  C  be  a  computational  structure  with  minimal  intricacy  realization  (A01  B01  CQ,  D0) 
with 


0>  = 


Cn 
C‘2 1 


C\2 

C22 


5 


Cf  1  invertible.  By  Lemma  3,  the  subsystem  structure  S  is  unique.  Since  Cu  is  invertible, 
we  see  that  equations  (58)  and  (60)  imply  that  the  minimal  intricacy  realization  uniquely 
specifies  the  dynamical  structure  function  of  the  system.  By  definition,  the  signal  structure 
is  unique.  Finally,  write  the  generalized  state-space  realization  of  C  as 


(A,B,C,D) 


A  A 
A  A 


B 

B 


,  [  C  C  ]  ,D 
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The  uniqueness  of  the  zero  pattern  structure  follows  from  its  one-to-one  correspondence  with 
the  transfer  function  G(s)  =  C0(sl  —  A0)_1i?0  +  Da  which  can  also  be  expressed  as 

(< C  +  C(I  -  A)~1A)(sI  -  A-A(I-  A)-1A)~1(kB  +  B(I  -  A)~lA). 


□ 

It  is  well  known  that  a  transfer  function  G(s)  can  be  realized  using  an  infinite  number 
of  state-space  realizations.  Without  additional  assumptions,  e.g.  full  state  feedback,  it  is 
impossible  to  uniquely  associate  a  single  state-space  realization  with  a  given  transfer  function. 
On  the  other  hand,  a  state  space  realization  specifies  a  unique  transfer  function.  In  this  sense, 
a  transfer  function  contains  less  information  than  the  state  space  realization. 

Similarly,  subsystem,  signal,  and  zero  pattern  structure  can  be  realized  using  multiple 
complete  computational  structures.  Without  additional  assumptions,  it  is  impossible  to 
associate  a  unique  complete  computational  structure  with  a  given  subsystem,  signal,  or  zero 
pattern  structure.  Theorem  1  shows  that  a  complete  computational  structure  specifies  a 
unique  subsystem,  signal,  and  zero  pattern  structure.  In  this  sense,  a  complete  computational 
structure  is  a  more  informative  description  of  system  structure  than  subsystem,  signal  and 
zero  pattern  structure.  The  next  result  has  a  similar  flavor  and  follows  directly  from  the 
one-to-one  correspondence  of  a  system’s  transfer  function  with  its  zero  pattern  structure. 

Theorem  2.  Every  subsystem  structure  or  signal  structure  specifies  a  unique  zero  pattern 
structure. 


Proof.  Consider  the  LFT  representation  E(N,  S )  of  a  subsystem  structure  S;  write 


N 


■  0 

I  ' 

L 

K 

as  in  equation  (54).  The  linear  fractional  transform  gives  the  input-output  map,  i.e.  the 
transfer  function.  Thus,  G(s)  =  (/  —  SK)~1SL. 

Similarly,  using  the  dynamical  structure  representation  of  the  signal  structure  W  given 
in  equation  (60),  we  can  solve  for  the  transfer  function 


G(s) 


(i- 

'  Q(s)  0  ' 

n 

\ 

C2  0 

>  i 

P(s)  +  (I-Q(s))D1 
D2 


The  result  follows  from  the  definition  of  zero  pattern  structure. 


□ 


The  relationship  between  subsystem  structure  and  signal  structure  is  not  so  straightfor¬ 
ward.  Nevertheless,  subsystem  structure  does  specify  a  unique  signal  structure  for  a  class  of 
systems,  namely  systems  with  subsystem  structure  composed  of  single  output  (SO)  subsys¬ 
tems  and  where  every  manifest  variable  is  involved  in  subsystem  interconnection.  For  this 
class  of  systems,  subsystem  structure  is  a  richer  description  of  system  structure  than  signal 
structure. 
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Single  Output  Subsystem  Structure  and  Signal  Structure 

Theorem  3.  Let  S  be  a  SO  subsystem  structure  with  LFT  representation  T (N ,  S) .  Suppose 


0 

I  ' 

L 

K 

where  [L  K~\  has  full  column  rank.  Then  S  uniquely  specifies  the  signal  structure  of  the 
system. 


Proof.  The  dynamics  of  N  and  S  satisfy 

V  = 

[0]  U  +  [/]  Y 

(68) 

7T  = 

LU  +  KY 

(69) 

V  = 

Sn 

(70) 

Combining  the  second  and  third  equation,  we  get 


Y  =  Sti  =  S[K  L ] 


Y 

U 


Since  S  is  a  SO  subsystem  structure,  the  entries  of  S  describe  direct  causal  dependencies 
among  manifest  variables  involved  in  interconnection.  Since  [L  iC]  has  full  column  rank 
and  is  a  binary  matrix,  this  means  that  each  manifest  variable  is  used  at  least  once  in 
interconnection.  Thus,  S  describes  the  direct  causal  dependencies  between  all  manifest 
variables  of  the  system  and  specifies  the  signal  structure  of  the  system.  □ 

Notice  that  for  the  class  of  systems  described  above,  the  four  representations  of  structure 
can  be  ordered  in  terms  of  information  content.  Theorem  1  shows  that  the  complete  com¬ 
putational  structure  uniquely  specifies  all  the  other  representations  of  structure  and  thus 
is  the  most  informative  of  the  four.  By  Theorem  3  and  Theorem  2  respectively,  subsystem 
structure  uniquely  specifies  the  signal  structure  and  zero  pattern  structure  of  the  system  and 
thus  is  the  second  most  informative.  Similarly,  signal  structure  is  the  third  most  informative 
and  zero  pattern  structure  is  the  least  informative  of  the  four  representations  of  structure. 

We  note  that  the  converse  of  Theorem  is  also  true,  namely  if  the  subsystem  structure  of 
a  system  specifies  a  unique  signal  structure  then  the  subsystem  structure  is  a  SO  subsystem 
structure  where  every  manifest  variable  is  an  interconnection  variable.  The  proof  is  simple 
and  follows  from  the  result  of  the  next  subsection.  We  also  provide  several  examples  that 
show  how  a  multiple  output  subsystem  structure  can  be  consistent  with  multiple  signal 
structures  -  these  all  will  serve  to  illustrate  the  general  relationship  between  subsystem  and 
signal  structure  outside  of  the  class  of  systems  mentioned  above. 


The  Relationship  Between  Subsystem  and  Signal  Structure  Subsystem  structure 
and  signal  structure  are  fundamentally  different  descriptions  of  system  structure.  In  general, 
subsystem  structure  does  not  encapsulate  the  information  contained  in  signal  structure.  Sig¬ 
nal  structure  describes  direct  causal  dependencies  between  manifest  variables  of  the  system. 
Subsystem  structure  describes  closed  loop  dependencies  between  manifest  variables  involved 
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in  the  interconnection  of  subsystems.  Both  representations  reveal  different  perspectives  of  a 
system’s  structure.  The  next  result  makes  this  relationship  between  subsystem  and  signal 
structure  precise. 


Theorem  4.  Given  a  system  G,  let  be  the  LFT  representation  of  a  subsystem 

structure  S.  In  addition,  let  the  signal  structure  of  the  system  G  be  denoted  as  in  equation 
(60).  Let  Y (Si)  denote  the  outputs  associated  with  subsystem  Si.  Define 


[Qint  (s).  ij 


Qij(s)  yi,yj  £  Y(Sk),Sk  a  subsystem  in  S 
0  otherwise, 


and  Qext  =  Q(s)  —  Qint(s).  Then  the  signal  structure  and  subsystem  structure  are  related  in 
the  following  way: 

S[  L  |  K  ]  =  (I -Qim)-1  [  P  |  Qext  ]  (71) 

Proof.  Examining  relation  (71),  observe  that  the  ijth  entry  of  the  left  hand  side  describes 
the  closed  loop  causal  dependency  from  the  jth  entry  of  [  UT  YT  ]T  to  Yt.  By  closed  loop, 
we  mean  that  they  do  not  describe  the  internal  dynamics  of  each  subsystem,  e.g.  the  direct 
causal  dependencies  among  outputs  of  a  single  subsystem.  Thus,  these  closed  loop  causal 
dependencies  are  obtained  by  solving  out  the  intermediate  direct  causal  relationships,  i.e. 
the  entries  in  Qint.  Notice  that  the  right  hand  side  of  (71)  also  describes  the  closed  loop  map 
from  [ UT  YT ] T  to  Y,  and  in  particular  the  i jth  entry  of 

(I  -  Qint)-1  [  P  |  Qext  } 


describes  the  closed  loop  causal  dependency  from  the  jth  entry  of 


'  U 


Y  f  to  Yt. 


□ 


As  a  special  case,  notice  that  for  SO  subsystem  structures,  Qint  becomes  the  zero  matrix 
and  that  for  subsystem  structures  with  a  single  subsystem,  S  becomes  the  system  transfer 
function,  L  becomes  the  identity  matrix,  Qint  =  Q,  and  Qext  and  K  are  both  zero  matrices. 
The  primary  import  of  this  result  is  that  a  single  subsystem  structure  can  be  consistent  with 
two  or  more  signal  structures  and  that  a  single  signal  structure  can  be  consistent  with  two 
or  more  subsystem  structures.  Consider  the  following  examples: 


Example  5.  A  Signal  Structure  consistent  with  two  Subsystem  Structures 
In  this  example,  we  will  show  how  a  signal  structure  can  be  consistent  with  two  subsystem 
structures.  To  do  this  we  construct  two  different  generalized  state-space  realizations  that  yield 
the  same  minimal  intricacy  realization  but  different  admissible  partitions  (see  Definition  7). 
The  result  is  two  different  subsystem  structures  that  are  associated  with  the  same  signal 
structure.  First,  we  consider  the  complete  computational  structure  C\  with  generalized  state- 
space  realization 


(Ai,  Bi,  Cl5  Dx) 


Ai  Ai 
Ai  Ai 


Bi 


,  [  Cj  Cx  ],D1 
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“j 


(a)  The  complete  computational  structure 
C  x  with  generalized  state-space  realization 
(Ai,Bi,  Ci,Di).  We  have  omitted  the  self-loops 
on  each  of  the  /  vertices  for  the  sake  of  visual 
clarity. 


(b)  The  complete  computational  structure 
C2  with  generalized  state-space  realization 
(A2,  B2,  C2,  D2).  We  have  omitted  the  self- loops 
on  each  of  the  /  vertices  for  the  sake  of  visual 
clarity. 


Figure  29:  The  complete  computational  structures  of  two  systems  that  differ  only  by  a  slight 
rearrangement  of  their  state-space  dynamics.  The  rearrangement  results  in  two  different 
subsystem  structure  representations.  Nonetheless,  the  resulting  minimal  intricacy  realization 
for  both  systems  is  the  same,  implying  that  both  systems  have  an  identical  signal  structure. 
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C\  —  04X5,  C 1  —  -Z4, 


and  D\  —  04Xi-  Figure  29a  shows  the  computational  structure  C\.  Next  consider  the  complete 
computational  structure  C2  with  generalized  state-space  realization 


(a2,  b2,  c2,  d2) 
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where 
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A2  —  A\,  A2  —  [0]4  , 


B2  —  B\  —  15x1, 


B2  —  Bi  —  04xi, 


C2  =  Cu 


c2  =  c1? 


and  D2  —  Di.  Figure  29b  shows  the  computational  structure  C2.  The  difference  between  these 
two  computational  structures  is  evident  more  in  the  subsystem  structure  representation  of 
the  system  -  note  how  replacing  Ai  with  A 2,  essentially  externalizes  internal  dynamics.  The 
result  is  that  C2  admits  a  subsystem  structure  S2  which  divides  one  of  the  subsystems  of  Si 
into  two  subsystems.  This  is  more  apparent  in  the  LFT  representations  of  Si  and  S2 ;  the 
LFT  representation  of  Si  is  given  by  fF(Ni,Si)  = 


and 


with  Sn,  S12,  S13  given  by 
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The  LFT  representation  of  S2  is  represented  as  the  LFT  J7(N2,S2)  =  where 
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However,  if  we  consider  the  minimal  intricacy  realizations  of  C\,C2  we  get  the  same 
state-space  realization  (A0,  Ba,  CQ,  Da )  with 
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The  signal  structure  of  the  system  is  thus  specified  by  the  dynamical  structure  function 
(■ Q,P)(s ),  with 
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Example  6.  A  Subsystem  Structure  consistent  with  two  Signal  Structures 
Now  we  consider  a  subsystem  structure  with  multiple  signal  structures.  Recalling  the  discus¬ 
sion  above,  subsystem  structure  describes  the  closed  loop  causal  dependencies  between  mani¬ 
fest  interconnection  variables  and  signal  structure  specifies  direct  causal  dependencies  between 
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manifest  variables.  When  it  is  impossible  to  determine  these  direct  causal  dependencies  by 
inspection  from  the  closed  loop  causal  dependencies  in  a  subsystem  structure  representation, 
then  there  can  be  multiple  signal  structures  that  are  consistent  with  the  subsystem  structure. 
Reconsider  S2-  The  LFT  is  given  by  P(N2,  S2)  where 
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with  S21  and  S22  given  by 
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respectively. 

If  we  consider  the  relation 


S2[K  L]  =  (/  -  QintY1  [  Qext  P], 
we  can  take  ( Q,P)(s )  to  equal 
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and  Q ext  —  Q 1  Qinti  or  Q2^s') 
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and  taking  Qint  =  [0]  and  Qext  =  Q  2-  It  is  routine  to  show  that  with  these  definitions,  both 
(Qi,Pi)(s)  and  ( Q2,P2)(s )  are  consistent  with  S2.  Thus  a  single  subsystem  structure  can  be 
consistent  with  two  signal  structures.  It  is  important  to  note  that  one  of  our  signal  structures 
exhibited  direct  causal  dependencies  corresponding  exactly  with  the  closed  loop  dependencies 
described  by  subsystem  structure.  We  note  that  this  correspondence  sometimes  occurs  (as  in 
the  case  of  Theorem  3)  but  generally  does  not  hold. 

4.5.4  Impact  on  Systems  Theory 

The  introduction  of  partial  structure  representations  suggests  new  problems  in  systems  the¬ 
ory.  These  problems  explore  the  relationships  between  different  representations  of  the  same 
system,  thereby  characterizing  properties  of  the  different  representations.  For  example,  clas¬ 
sical  realization  theory  considers  the  situation  where  a  system  is  specified  by  a  given  transfer 
function,  and  it  explores  how  to  construct  a  consistent  state  space  description.  Many  im¬ 
portant  ideas  emerge  from  the  analysis: 

1.  State  realizations  are  generally  more  informative  than  a  transfer  function  representa¬ 
tion  of  a  system,  as  there  are  typically  many  state  realizations  consistent  with  the  same 
transfer  function, 

2.  Order  of  the  state  realization  is  a  sensible  measure  of  complexity  of  the  state  represen¬ 
tation,  and  there  is  a  well-defined  minimal  order  of  any  realization  consistent  with  a 
given  transfer  function;  this  minimal  order  is  equal  to  the  Smith-McMillian  degree  of 
the  transfer  function, 

3.  Ideas  of  controllability  and  observability  of  a  state  realization  characterize  important 
properties  of  the  realization,  and  any  minimal  realization  is  both  controllable  and 
observable. 

In  a  similar  way,  introducing  partial  structure  representations  impacts  a  variety  of  concepts 
in  systems  theory,  including  realization,  minimality,  and  model  reduction. 
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Realization  The  definition  of  partial  structure  representations  enrich  the  kinds  of  real¬ 
ization  questions  one  may  consider.  In  the  previous  section,  we  demonstrated  that  partial 
structure  representations  of  a  system  are  generally  more  informative  than  its  transfer  function 
but  less  informative  than  its  state  realization.  Thus,  two  classes  of  representation  questions 
emerge:  reconstruction  problems  and  structure  realization  problems  (Figure  30). 


Figure  30:  Partial  structure  representations  introduce  new  classes  of  realization  problems: 
reconstruction  and  structure  realization.  These  problems  are  distinct  from  identification, 
and  each  depends  on  the  type  of  partial  structure  representation  considered. 

Reconstruction  problems  consider  the  construction  of  a  partial  structure  representation  of 
a  system  given  its  transfer  function.  Because  partial  structure  representations  are  generally 
more  informative  than  a  transfer  function,  these  problems  are  ill-posed-like  the  classical 
realization  problem.  Nevertheless,  one  may  consider  algorithms  for  generating  canonical 
representations,  or  one  may  characterize  the  additional  information  about  a  system-beyond 
knowledge  of  its  transfer  function-that  one  would  require  in  order  to  specify  one  of  its 
partial  structure  representations.  In  particular,  we  may  consider  two  types  of  reconstruction 
problems: 

1.  Signal  Structure  Reconstruction:  Given  a  transfer  function  G(s)  with  its  associ¬ 
ated  zero-pattern  structure  Z ,  find  a  signal  structure,  W,  with  dynamical  structure 
function  (Q,  P )  as  in  Equation  (60)  such  that  G  =  (/  —  Q)~lP , 

2.  Subsystem  Structure  Reconstruction:  Given  a  transfer  function  G(s)  with  its 
associated  zero-pattern  structure  Z,  find  a  subsystem  structure,  S,  with  structured 
LFT  (TV,  S )  as  in  Equation  (54)  such  that  G  =  (/  —  SK )~1SL. 

Signal  structure  reconstruction  is  also  called  network  reconstruction,  particularly  in  systems 
biology  where  it  plays  a  central  role.  There,  the  objective  is  to  measure  fluctuations  of  various 
proteins,  or  other  chemical  species,  in  response  to  particular  perturbations  of  a  biochemical 
system,  and  then  infer  causal  dependencies  among  these  species  [78,  79,  76,  71,  80,  81]. 

Structure  realization  problems  then  consider  the  construction  of  a  state  space  model, 
possibly  generalized  to  include  auxiliary  variables  as  necessary,  consistent  with  a  given  partial 
structure  representation  of  a  system.  Like  the  classical  realization  problem  or  reconstruction 
problems,  these  problems  are  also  ill-posed  since  there  are  typically  many  state  realizations 
of  a  given  partial  structure  representation  of  a  system.  Also,  like  reconstruction  problems, 
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these  realization  problems  can  be  categorized  in  two  distinct  classes,  depending  on  the  type 
of  partial  structure  representation  that  is  given: 

3.  Signal  Structure  Realization:  Given  a  system  G  with  signal  structure  W  and 
associated  dynamical  structure  function  (Q,P),  find  a  state  space  model  (A,  B,C,  D) 
consistent  with  ( Q,P ),  i.e.  such  that  Equations  (57-60)  hold, 

4.  Subsystem  Structure  Realization:  Given  a  system  G  with  subsystem  structure  S 
and  associated  structured  ( N,S )  (recorded  in  the  LFT  representation  T(N,  S)),  find 
a  generalized  state  space  model  of  the  form  of  Equation  (48)  consistent  with  (. N ,  S). 

Signal  structure  realization  may  sometimes  be  called  network  realization ,  consistent  with  the 
nomenclature  for  signal  structure  reconstruction. 

Note  that  all  the  reconstruction  and  structure  realization  problems  here  are  distinct  from 
identification  problems,  just  as  classical  realization  is  different  from  identification.  For  the 
systems  considered  here,  identification  refers  to  the  use  of  input-output  data  (and  no  other 
information  about  a  system)  to  choose  a  representation  that  best  describes  the  data  in  some 
appropriate  sense.  Because  input-output  data  only  characterizes  the  input-output  map  of  a 
system,  identification  can  at  best  characterize  the  system’s  transfer  function;  no  information 
about  structure,  beyond  the  zero-pattern  of  the  transfer  function,  is  available  in  such  data. 
In  spite  of  this  distinction,  however,  it  is  not  uncommon  for  reconstruction  problems  to  be 
called  structure  identification  problems.  Nevertheless,  one  may  expose  such  problems  as  the 
concatenation  of  an  identification  and  a  reconstruction  problem  and  precisely  characterize 
the  extra  information  needed  to  identify  such  structure  by  carefully  analyzing  the  relevant 
reconstruction  problem,  independent  of  the  particular  identification  technique  [71,  78]. 

Minimality  Just  as  partial  structure  representations  enrich  the  classical  realization  prob¬ 
lem,  they  also  impact  the  way  we  think  about  minimality.  Certainly  the  idea  of  a  minimal 
complexity  representation  is  relevant  for  each  of  the  four  problems  listed  above,  but  clearly 
the  relevant  notion  of  complexity  may  be  different  depending  on  the  representation.  We 
consider  each  situation  as  follows: 

1.  Minimal  Signal  Structure  Reconstruction:  In  this  situation  one  needs  to  consider 
how  to  measure  the  complexity  of  a  system’s  signal  structure,  W,  or  its  associated 
dynamical  structure  function,  (Q,P).  Some  choices  one  may  consider  could  be  the 
number  of  edges  in  W,  the  maximal  order  of  any  element  in  ( Q,P ),  or  the  maximal 
order  of  any  path  from  any  input  Ui  to  any  output  yj.  The  problem  would  then  be  to 
find  a  minimal  complexity  signal  structure  consistent  with  a  given  transfer  function. 

2.  Minimal  Subsystem  Structure  Reconstruction:  In  this  situation  one  needs  to 
consider  how  to  measure  the  complexity  of  a  system’s  subsystem  structure,  S,  or  its 
associated  structured  LFT,  ( N ,  S).  One  notion  could  be  to  measure  complexity  by  the 
number  of  distinct  subsystems;  the  problem  would  then  be  to  find  the  minimal  com¬ 
plexity  subsystem  structure  representation  consistent  with  a  given  transfer  function. 
Another  notion  could  be  the  number  of  non-zero  entries  in  the  S  matrix,  where  (N,  S ) 
denote  the  the  LFT  associated  with  the  subsystem  structure.  Using  this  measure, 
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a  single  subsystem  block  with  no  zero  entries  would  be  considered  a  more  complex 
representation  than  a  subsystem  structure  with  a  large  number  of  distinct,  albeit  in¬ 
terconnected,  subsystems. 

3.  Minimal  Signal  Structure  Realization:  In  this  situation  one  needs  to  consider 
how  to  measure  the  complexity  of  a  system’s  zero-intricacy  state  realization,  from 
which  signal  structure  is  derived.  The  obvious  choice  would  be  to  use  the  order  of  the 
realization  as  a  natural  measure  of  complexity,  and  the  problem  would  then  become 
to  find  the  minimal  order  state  realization  consistent  with  a  given  signal  structure,  W, 
or,  equivalently,  with  its  associated  dynamical  structure  function,  ( Q,P ).  Note  that 
this  minimal  order  is  guaranteed  to  be  finite  (for  the  systems  considered  here)  and  can 
easily  be  shown  to  be  greater  or  equal  to  the  Smith-McMillian  degree  of  the  transfer 
function  specified  by  the  signal  structure;  we  call  this  number  the  structural  degree  of 
the  signal  structure  [72], 

4.  Minimal  Subsystem  Structure  Realization:  In  this  situation  one  needs  to  con¬ 
sider  how  to  measure  the  complexity  of  a  generalized  state  realization  in  the  presence 
of  auxiliary  variables.  Both  the  order  and  the  intricacy  of  the  realization  offer  some 
perspective  of  its  complexity,  but  one  needs  to  consider  how  they  each  contribute  to 
the  overall  complexity  of  the  realization.  The  problem  would  then  be  to  find  a  minimal 
complexity  generalized  state  realization  consistent  with  a  given  subsystem  structure. 

These  various  problems  demand  new  ideas  for  thinking  about  the  complexity  of  a  sys¬ 
tem’s  representation,  especially  that  of  a  partial  structure  representation.  These  new  ideas 
about  complexity,  in  turn,  introduce  opportunities  for  characterizing  minimality  of  a  rep¬ 
resentation  in  terms  that  add  insight  to  our  understanding  of  the  relationship  between  a 
system’s  behavior  and  its  structure,  much  as  controllability  and  observability  characterize 
classical  notions  of  minimality  in  a  system’s  state  realization.  Besides  suggesting  the  need  for 
a  characterization  of  minimality,  however,  these  ideas  also  impact  notions  of  approximation 
and  how  we  think  about  model  reduction. 

Model  Reduction  Each  of  the  reconstruction  and  structure  realization  problems  de¬ 
scribed  above  have  associated  with  them  not  only  a  minimal-representation  version  of  the 
problem,  but  also  an  approximate-representation  version  of  the  problem.  The  minimal- 
representation  versions  of  these  problems,  as  described  above,  seek  to  construct  a  represen¬ 
tation  of  minimal  complexity  in  the  targeted  class  that  is  nevertheless  consistent  with  the 
system  description  provided.  Similarly,  approximate-representation  versions  of  these  prob¬ 
lems  seek  to  construct  a  representation  in  the  targeted  class  that  has  a  lower  complexity 
than  the  minimal  complexity  necessary  to  deliver  consistency  with  the  system  description 
provided.  As  a  result,  consistency  with  the  given  system  description  can  not  be  achieved,  so 
measures  of  approximation  become  necessary  to  sensibly  discuss  a  “best”  representation  of 
the  specified  complexity. 

For  example,  associated  with  the  classical  realization  problem  is  the  standard  model 
reduction  problem.  In  this  situation,  a  transfer  function  is  specified,  and  one  would  like  to 
construct  a  state  realization  with  a  complexity  that  is  lower  than  that  which  is  necessary 
(for  such  a  realization  to  be  consistent  with  the  given  transfer  function)  that  nevertheless 
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“best”  approximates  it.  Note  that  the  appropriate  notion  of  complexity  depends  on  the  target 
representation  class;  here,  the  target  representation  class  is  the  set  of  state  realizations,  so  the 
appropriate  notion  of  complexity  is  model  order.  Likewise,  note  that  the  appropriate  notion 
of  approximation  depends  on  the  type  of  system  representation  that  is  initially  provided; 
here,  a  transfer  function  is  provided,  so  an  appropriate  measure  of  approximation  could  be 
an  induced  norm,  such  as  H^.  Thus,  one  could  measure  the  quality  of  an  approximation 
by  measuring  the  induced  norm  of  the  error  between  the  given  transfer  function  and  that 
specified  by  the  approximate  state  realization.  Note  that  since  the  H <*,  model  reduction 
problem  remains  open,  many  alternative  measures  and  approaches  have  been  suggested. 
In  any  event,  because  the  specified  system  description  is  a  transfer  function,  the  resulting 
measure  of  approximation  is  typically  one  that  either  directly  or  indirectly  measures  the 
difference  in  input-output  dynamic  behavior  between  the  approximate  model  and  the  given 
system  description;  the  focus  is  on  dynamics,  not  system  structure,  when  considering  notions 
of  approximation  in  the  standard  model  reduction  problem. 

Similarly,  there  is  an  appropriate  reduction  problem  associated  with  each  of  the  minimal 
reconstruction  and  minimal  structure  realization  problems  described  above.  Like  standard 
model  reduction,  the  appropriate  notion  of  complexity  is  characterized  by  the  class  of  target 
representations,  while  the  appropriate  measure  of  approximation  depends  on  the  system 
representation  that  is  initially  provided.  To  make  these  ideas  more  concrete,  we  discuss  each 
of  the  problems  individually: 

1.  Approximate  Signal  Structure  Reconstruction:  In  this  situation  one  would  like 
to  find  a  signal  structure  representation  with  lower  complexity  e.g.  fewer  number  of 
edges,  etc.  than  the  minimal  level  of  complexity  necessary  to  specify  a  given  transfer 
function.  When  using  the  number  of  edges  as  a  complexity  measure,  this  problem 
may  be  interpreted  as  trying  to  restrict  the  number  of  communication  links  between 
individual  computation  nodes  in  a  distributed  system  while  approximating  some  de¬ 
sired  global  dynamic.  The  appropriate  measure  of  approximation,  then,  is  any  measure 
that  sensibly  compares  the  desired  transfer  function  from  that  specified  by  the  reduced 
signal  structure,  such  as  H0 Q. 

2.  Approximate  Subsystem  Structure  Reconstruction:  In  this  situation  one  would 
like  to  find  a  subsystem  structure  representation  with  lower  complexity  than  the  mini¬ 
mal  level  of  complexity  necessary  to  specify  a  given  transfer  function.  If  one  considers 
the  number  of  subsystems  as  the  measure  of  complexity,  then  this  problem  is  trivial 
since  any  transfer  function  is  consistent  with  a  subsystem  structure  with  a  single  sub¬ 
system.  If  one  measures  complexity  by  the  number  of  non-zero  entries  in  the  S  matrix, 
where  (A,  S)  denote  the  the  LFT  associated  with  the  subsystem  structure,  then  it 
appears  likely  that  one  could  often  formulate  a  meaningful  approximation  problem. 

Unlike  these  approximate  reconstruction  problems,  approximate  realization  problems 
may  have  dual  relevant  measures  of  approximation,  since  having  an  initial  partial  struc¬ 
ture  representation  of  a  system  also  always  specifies  its  transfer  function.  Measuring  the 
similarity  between  transfer  functions  determines  the  degree  to  which  a  lower  order  system 
approximates  the  dynamics  of  a  given  system,  while  measuring  the  similarity  between  partial 
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structure  representations  determines  the  degree  to  which  a  lower  order  system  approximates 
the  structure  of  a  given  system. 

3.  Approximate  Signal  Structure  Realization:  In  this  situation  one  would  like  to 

find  a  state  realization  with  a  model  complexity  that  is  lower  than  the  minimal  com¬ 
plexity  necessary  to  specify  a  given  signal  structure.  Typically  the  complexity  of  a 
generalized  state  realization  would  be  measured  by  both  the  order  and  intricacy  of  the 
realization.  Since  signal  structure  only  depends  on  the  zero  intricacy  realization  of  a 
system,  however,  a  minimal  complexity  realization  will  always  have  zero  intricacy;  thus 
model  order  is  the  only  relevant  notion  of  complexity.  Moreover,  since  the  structural 
degree  of  a  particular  signal  structure  may  be  strictly  greater  than  the  Smith-McMillian 
degree  of  the  transfer  function  it  specifies,  a  few  distinct  kinds  of  approximation  prob¬ 
lems  emerge:  structural  approximation  and  dual  approximation  (Figure  31). 

(a)  Structural  Approximation:  When  the  order  of  the  approximation  is  less  than 
the  structural  degree  of  the  given  signal  structure  but  greater  than  the  Smith- 
McMillian  degree  of  the  associated  transfer  function,  the  resulting  approximation 
can  specify  the  given  transfer  function  exactly,  even  though  its  signal  structure 
only  approximates  that  of  the  given  system.  Note  that  the  sense  in  which  simi¬ 
larity  in  signal  structure  should  be  measured  is  an  area  of  on-going  research. 

(b)  Dual  Approximation:  When  the  order  of  the  approximation  is  less  than  the 
Smith-McMillian  degree  of  the  transfer  function  specified  by  the  given  signal  struc¬ 
ture,  then  it  will  represent  neither  the  structure  nor  the  dynamics  of  the  given 
system  exactly. 

4.  Approximate  Subsystem  Structure  Realization:  In  this  situation  one  would  like 

to  find  a  generalized  state  realization  with  a  model  complexity  that  is  lower  than  the 
minimal  complexity  necessary  to  specify  a  given  subsystem  structure.  Here,  complex¬ 
ity  of  a  generalized  state  realization  is  measured  both  in  terms  of  intricacy  and  order 
since  intricacy  of  a  realization  directly  impacts  the  number  of  admissible  subsystem 
blocks  (see  Figure  23)  while  order  impacts  the  ability  of  the  realization  to  approximate 
the  transfer  function  specified  by  the  given  subsystem  structure.  As  a  result,  three  dis¬ 
tinct  approximation  problems  emerge  to  complement  subsystem  structure  realization: 
Structure-Preserving  Model  Reduction,  Subsystem  Structure  Approximation,  Subsys¬ 
tem  Dual  Approximation  (Figure  32). 

(a)  Structure-Preserving  Model  Reduction:  When  the  intricacy  of  an  approx¬ 
imation  is  high  enough,  the  subsystem  structure  of  a  system  can  be  preserved 
while  its  dynamic  behavior  is  approximated  by  lower-order  systems.  A  naive  ap¬ 
proach  to  such  reduction  would  be  to  reduce  the  order  of  various  subsystems. 
However,  even  when  each  subsystem  is  well  approximated,  the  dynamic  behavior 
of  the  overall  interconnection  may  be  very  different  from  the  original  system.  As  a 
result,  a  rich  literature  has  developed  in  this  area  that  develops  sophisticated  tech¬ 
niques  for  approximating  the  dynamics  of  the  closed-loop,  interconnected  system 
while  preserving  its  subsystem  structure  [82,  83]. 
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Figure  31:  Approximate  signal  structure  realization  leads  to  two  distinct  types  of  reduction 
problems.  Structural  approximation  exactly  realizes  the  dynamics  of  a  given  signal  structure 
while  only  approximating  its  structure;  dual  approximation  captures  neither  the  dynamics 
nor  the  structure  of  a  given  signal  structure  exactly. 

(b)  Subsystem  Structure  Approximation:  When  the  order  of  an  approximation 
is  high  enough,  the  dynamic  behavior  of  a  system  can  be  preserved  while  is  subsys¬ 
tem  structure  is  approximated  by  lower-intricacy  realizations.  The  sense  in  which 
similarity  of  subsystem  structure  should  be  measured  remains  an  open  topic  for 
research. 

(c)  Subsystem  Dual  Approximation:  When  both  the  order  and  the  intricacy  of 
an  approximation  are  lower  than  the  minimal  values  necessary  to  realize  a  given 
subsystem  structure  and  the  transfer  function  it  specifies,  then  the  objective  of 
the  reduction  problem  is  to  find  the  realization  of  the  specified  complexity  that 
best  approximates  the  structure  and  dynamics  of  the  given  system. 

The  introduction  of  partial  structure  representations  suggests  a  number  of  new  prob¬ 
lems  in  systems  theory.  These  problems  include  new  classes  of  realization  problems,  called 
reconstruction  and  structural  realization  problems,  as  well  as  a  number  of  new  reduction 
problems.  Each  of  these  problems  differ  depending  on  the  partial  structure  representation 
one  considers,  and  a  number  of  research  issues  remain  to  properly  formulate  most  of  them. 
The  overview  offered  here  is  merely  meant  to  give  a  perspective  of  the  landscape  of  problems 
that  emerges  with  the  introduction  of  partial  structure  representations. 
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Figure  32:  Approximate  subsystem  structure  realization  leads  to  three  distinct  types  of 
reduction  problems:  1)  Structure-preserving  model  reduction  preserves  the  subsystem  struc¬ 
ture  of  a  given  system  while  approximating  its  dynamic  behavior,  2)  Subsystem  structure 
approximation  preserves  the  dynamic  behavior  but  approximates  the  subsystem  structure, 
and  3)  Subsystem  dual  approximation  captures  neither  the  structure  nor  the  dynamic  be¬ 
havior  of  a  system  described  by  a  particular  subsystem  structure. 
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4.5.5  Summary 

This  section  introduced  the  idea  of  a  system’s  complete  computational  structure  as  a  base¬ 
line  for  comparing  simplified  representations.  Although  closely  aligned  with  a  system’s  state 
space  realization,  we  demonstrated  the  need  for  auxiliary  variables  to  encode  information 
about  a  system’s  admissible  subsystem  structure.  This  results  in  a  “generalized”  state  de¬ 
scription  that  is  a  differential  algebraic  system  with  differentiation  index  zero,  and  the  num¬ 
ber  of  auxiliary  variables  becomes  an  additional  measure  of  complexity  (besides  model  order) 
that  we  call  intricacy.  This  generalized  state  representation,  and  its  associated  graphical  rep¬ 
resentation,  is  the  system’s  complete  computational  structure. 

Partial  structure  representations  were  then  introduced  as  a  means  for  simplifying  a  sys¬ 
tem’s  structural  description  while  retaining  a  complete  representation  of  its  input-output  dy¬ 
namic  behavior.  These  included  subsystem  structure,  signal  structure,  and  the  zero-pattern 
of  the  system’s  transfer  function. 

Subsystem  structure  was  first  introduced  as  the  most  refined  view  of  the  interconnection 
of  a  system’s  legitimate  subsystems.  This  description  is  represented  as  the  lower  linear 
fractional  transformation  of  a  static  interconnection  matrix  N  with  a  block  diagonal  operator 
S.  It’s  graphical  representation  is  a  block  diagram,  like  the  system’s  complete  computational 
structure,  that  is  a  condensation  graph  with  respect  to  a  meaningful  partition  of  the  system’s 
states. 

Signal  structure,  on  the  other  hand,  is  a  signal  flow  diagram  of  the  causal  dependen¬ 
cies  among  manifest  variables  given  by  the  dynamical  structure  function  of  the  system.  We 
demonstrated  that  systems  may  exhibit  extremely  structured  behavior  in  the  signal  struc¬ 
ture  sense  while  having  no  apparent  structure  in  the  subsystem  or  computational  structure 
sense.  Moreover,  the  transformations  between  manifest  variables,  represented  as  edges  in  the 
signal  structure  graph,  do  not  necessarily  partition  the  system  states,  as  do  the  nodes  of  the 
subsystem  structure.  This  fact  implies  that  the  minimal  order  to  realize  a  particular  signal 
structure  may  be,  in  fact,  higher  than  the  minimal  order  necessary  to  realize  the  transfer 
function  specified  by  the  given  signal  structure. 

Finally,  the  weakest  notion  of  structure  is  the  pattern  of  zero  entries  in  the  system’s  trans¬ 
fer  function  matrix,  which  graphically  also  becomes  a  signal  flow  graph  like  signal  structure. 
We  demonstrated  that  this  representation  is  very  weak,  reminding  readers  that  a  diagonal 
transfer  function  matrix,  for  example,  does  not  imply  that  even  a  minimal  realization  of  the 
system  is  necessarily  decoupled.  Thus,  this  notion  of  structure  really  only  makes  a  statement 
about  the  closed-loop  dependencies  of  inputs  on  outputs  of  the  system. 

These  representations  were  then  shown  to  contain  differing  levels  of  structural  infor¬ 
mation  about  a  system.  In  particular,  it  was  shown  that  the  complete  computational  struc¬ 
ture  uniquely  specifies  both  the  subsystem  and  signal  structure  of  a  system,  and  that  either 
of  these  partial  structure  representations  uniquely  specify  the  transfer  function  (and  thus 
its  associated  zero-pattern  structure)  of  the  system.  Nevertheless,  the  relationship  between 
subsystem  and  signal  structure  is  less  definitive,  as  we  demonstrated  that  two  realizations 
of  the  same  system  may  share  subsystem  structure  but  have  different  signal  structures,  or, 
conversely,  two  realization  of  the  same  system  may  share  signal  structure  but  have  differ¬ 
ent  subsystem  structures.  These  different  representations  simply  appear  to  offer  different 
perspectives  of  the  system’s  structural  properties. 
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We  then  surveyed  the  landscape  of  new  problems  in  systems  theory  that  arise  when  con¬ 
sidering  these  partial  structure  representations.  In  particular,  we  showed  that  the  standard 
realization  problem  now  becomes  two  types  of  problems,  a  reconstruction  problem,  where  a 
transfer  function  is  given  and  one  would  like  to  determine  a  partial  structure  representation 
that  is  consistent  with  it,  and  a  structure  realization  problem,  where  a  partial  structure 
representation  is  given  and  one  would  like  to  find  a  generalized  state  realizations  that  is 
consistent  with  it.  Minimal  versions  of  these  problems  are  obtained  if  one  can  define  a 
sensible  notion  of  complexity  of  each  kind  of  representation,  and  various  suggestions  were 
offered.  Associated  approximation,  or  model  reduction  problems  were  then  characterized, 
where  the  target  representation  is  simpler  than  the  minimal  complexity  necessary  to  yield  a 
representation  that  is  consistent  with  the  given  model  or  description  of  the  system.  Here  we 
note  that  the  model  reduction  problems  begin  to  consider  structure  approximation  as  well  as 
approximation  of  the  system’s  dynamic  behavior,  leading  to  a  variety  of  new  problems  one 
may  consider.  It  is  out  hope  that  many  of  these  problems  will  be  addressed  in  the  coming 
years,  leading  to  a  more  thorough  and  complete  understanding  of  the  meaning  of  structure 
and  its  relationship  to  the  behavior  of  interconnected  dynamic  systems. 

4.6  Vulnerable  Links  and  Secure  Architectures  in  the  Stabilization 
of  Networks  of  Controlled  Dynamical  Systems 

The  Stuxnet  virus  attacked  an  Iranian  nuclear  power  plant  in  2010  and  caused  the  centrifuge’s 
rotors  to  malfunction  [98] .  It  gained  much  news  coverage  as  the  first  virus  attack  on  industrial 
systems.  Although  no  serious  damage  was  done,  it  has  highlighted  the  necessity  of  improving 
our  understanding  of  the  security  of  control  systems. 

Researchers  have  predicted  that  attacks  on  industrial  control  systems  would  increase 
[86,  84],  As  the  systems  are  becoming  more  networked,  securing  them  has  become  much 
harder  and  attacking  them  has  gotten  easier.  In  the  past,  securing  the  physical  plants  was 
enough  to  secure  the  system,  but  now  in  the  networked  architecture  the  communication 
channels  have  to  be  secured  too.  This  is  almost  impossible  to  achieve,  especially  when 
these  systems  are  being  connected  to  the  Internet,  with  connection  features  as  powerful 
as  remote  access  to  the  control  centers.  Regardless  of  the  improvement  in  the  industry’s 
secure  communications,  cryptography,  etc.,  a  simple  human  error  like  someone  forgetting  to 
change  a  default  password  on  their  account  could  give  an  attacker  complete  access  to  the 
resources  necessary  to  carry  out  a  complicated  attack.  Although  the  security  risks  are  real, 
these  systems  being  less  networked  in  the  future  is  highly  unlikely  because  of  the  usability 
advantages  that  a  networked  setting  offers. 

As  a  result,  considering  the  security  of  networked  control  systems  has  become  very  im¬ 
portant.  A  good  design  should  make  detecting  attacks  easy,  help  understand  the  effects  of 
an  attack,  make  it  difficult  to  execute  an  attack,  and  finally  minimize  the  consequences  if  an 
attack  is  successful.  This  section  will  contribute  in  designing  more  secure  networked  systems 
by  identifying  conditions  when  a  link  is  completely  secure  against  attacks  that  attempt  to 
destabilize  the  system.  Our  result  also  gives  a  measure  of  link  vulnerability,  which  corre¬ 
sponds  to  the  minimum  size  of  a  destabilizing  attack  on  the  link.  This  can  be  a  useful  tool 
to  understand  the  security  of  a  networked  system. 
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In  this  section,  we  view  security  as  a  robustness  issue  and  focus  on  making  systems  robust 
against  perturbations  on  a  single  communication  channel.  First,  we  give  a  summary  of  the 
types  of  attacks  that  a  system  might  suffer.  Then,  in  Section  7.6.4  we  present  our  main 
result.  Finally  in  Section  7.6.5  we  give  some  examples  to  illustrate  the  applications  of  our 
theory. 

4.6.1  Attack  Models 

In  the  literature,  attacks  on  control  systems  have  been  classified  into  two  types:  denial 
of  service  attacks,  when  the  attacker  jams  a  channel  in  order  to  destabilize  the  system, 
and  deception  attacks,  when  the  attack  adds  perturbations  on  particular  links  in  order  to 
compromises  the  reliability  of  the  controller’s  state  estimates  [96,  85].  We  consider  a  hybrid 
attack  model  where  the  attacker  adds  perturbations  to  the  channels,  not  just  jam  them,  in 
order  to  destabilize  the  system.  We  call  this  type  of  attack  a  destabilizing  attack. 

Denial  of  Service  (DoS)  Attack  Denial  of  service  attacks  prevent  signals  from  reaching 
their  intended  destination.  This  is  probably  the  easiest  and  most  common  attack,  and  it  is 
modeled  as  removal  of  an  edge  in  an  interconnected  structure.  It  might  be  done  by  jam¬ 
ming  the  communication  channel,  disrupting  the  transmitter/receiver,  changing  the  routing 
protocol,  saturating  the  receiver  with  extraneous  signals,  etc.  The  attacker’s  intent  of  such 
an  attack  could  be  to  degrade  the  system  performance  or  to  completely  destabilize  it.  [94] 
shows  that  performance  of  networked  control  systems  could  decrease  significantly  under  a 
DoS  attack.  [96]  gives  a  method  to  find  an  optimal  controller  that  minimizes  the  effect  of 
such  an  attack  on  linear  control  systems. 

In  [97],  the  authors  study  whether  a  DoS  attack  on  certain  links  can  make  the  system 
unreachable  or  uncontrollable.  They  also  develop  graph  theoretic  algorithms  to  identify  the 
minimal  number  of  edges  which  are  necessary  for  preserving  controllability  and  observability. 

Deception  Attack  The  goal  of  a  deception  attack  is  to  change  the  state  estimates  com¬ 
puted  by  a  model-based  controller.  This  type  of  attack  is  modeled  as  a  stable  additive 
perturbation  to  an  edge  in  the  network.  All  stabilizing  controllers  make  the  closed  loop 
system  stable,  hence,  a  stabilizing  controller  is  necessarily  stabilizable  from  the  plant.  So,  if 
at  attacker  gains  access  to  the  communication  channel  between  the  plant  and  the  controller, 
state  estimates  of  a  model-based  controller  can  be  altered.  To  prevent  this,  many  real  sys¬ 
tems  such  as  power  systems,  sensor  networks,  etc.,  are  equipped  with  a  Bad  Data  Detector 
(BDD)  [99,  92,  100].  A  BDD,  using  the  model  of  the  plant,  detects  deviation  of  the  state  es¬ 
timates  from  the  expected  and  raises  an  alarm  to  notify  the  human  operator.  Because  of  the 
presence  of  measurement  noise,  this  deviation  is  never  zero,  so  the  BDD  ignores  deviations 
that  are  smaller  than  a  specified  threshold.  Hence,  in  the  presence  of  BDDs,  the  attack  has 
to  change  the  state  estimates  without  increasing  the  chance  of  raising  an  alarm. 

In  [99]  the  authors  study  this  kind  of  attack  in  the  context  of  a  power  system.  They 
shown  that  it  is  in  fact  possible  for  an  attacker  to  change  the  state  estimates  to  a  specific 
value  without  increasing  the  chance  of  being  detected.  [100]  studies  a  similar  problem  in  the 
scenario  of  a  wireless  sensor  networks.  It  maps  an  approximation  of  the  set  of  all  possible 
values  the  attacker  could  drive  the  estimates  to. 
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[92]  studies  a  slightly  different  problem.  Here,  the  goal  of  the  attacker  is  to  change  the 
estimate  of  one  of  the  states  without  increasing  the  chance  of  being  detected.  The  authors 
recognize  that  while  doing  this  the  attacker  might  want  to  use  the  fewest  channels  possible 
or  might  try  to  keep  the  magnitude  of  the  attack  signal  small.  For  each  type  of  attack,  the 
authors  then  give  a  formulation  of  a  security  index  of  the  system. 

Destabilizing  Attack  Like  deception  attacks,  these  attacks  effectively  arise  as  an  additive 
perturbation  on  a  link  in  the  system  interconnection  structure.  Llnlike  deception  attacks, 
however,  they  seek  to  destabilize  the  system  rather  than  simply  move  the  system  state  to 
a  desired  value  without  being  detected.  BDDs  are  clearly  capable  of  detecting  the  destabi¬ 
lization  resulting  from  such  attacks,  nevertheless  serious  damage  and  even  complete  plant 
shut-down  may  result  by  the  time  system  operators  are  able  to  do  anything  about  it. 

A  rich  literature  in  systems  and  control  theory  explores  the  destabilization  of  systems  due 
to  additive  perturbations,  see  for  example  [91]  and  the  references  therein.  Security  analysis 
of  destabilizing  attacks  thus  appears  to  be  a  robustness  problem  with  respect  to  certain 
classes  of  perturbations.  Indeed,  we  adopt  this  point  of  view,  and  consider  security  problems 
to  be  essentially  robustness  problems  of  various  types. 

The  contribution  of  this  work,  then,  applied  to  this  class  of  attacks,  is  in  the  solution 
of  a  certain  class  of  robustness  problems  over  a  particular  kind  of  link  model-corresponding 
to  logical,  rather  than  the  physical,  links  of  a  system-and  with  respect  to  a  specific  class 
of  perturbations.  Llnlike  standard  robustness  measures  that  generally  consider  destabilizing 
perturbations  acting  over  all  channels  and  nodes  of  a  system,  here  we  restrict  our  attention 
specifically  to  perturbations  that  disrupt  a  single  link  in  the  system’s  signal  structure.  Our 
analysis  then  considers  such  single-link  perturbations  over  all  possible  system  links.  In  the 
next  section  we  explore  our  link  model  in  detail. 

4.6.2  Link  Models 

The  destabilizing  attacks  considered  here  are  additive  perturbations  acting  on  a  single  link  in 
a  system’s  logical  interconnection  structure.  There  are  many  characterizations  of  a  system’s 
structure,  see  for  example  [88,  89,  87].  One  characterization  would  consider  the  intercon¬ 
nection  structure  among  subsystems.  This  definition  of  structure,  also  called  the  system’s 
subsystem  structure,  would  represent  the  physical  interconnection  between  physical  compo¬ 
nents  of  a  particular  networked  system.  Linder  this  notion  of  structure,  a  link  would  represent 
the  signal  passing  between  two  subsystem  nodes  within  the  subsystem  interconnection  ar¬ 
chitecture.  In  contrast  to  the  subsystem  structure,  this  work  considers  another  definition  of 
system  structure  and,  consequently,  a  different  notion  of  a  system  link. 

In  this  work,  we  consider  a  partition  on  signals  of  the  system  into  two  categories:  exposed 
signals  and  hidden  signals.  The  logical  interconnection  structure,  or  architecture-also  called 
the  system’s  signal  structure-is  the  causal  relationship  between  exposed  signals  in  the  system. 
In  this  definition  of  structure,  a  link  is  a  system  describing  the  causal  dependency  between 
two  exposed  signal  nodes  of  the  logical  interconnection  architecture. 

Some  important  consequences  of  this  definition  of  link  include  the  fact  that  a  link  may 
represent  a  very  indirect  and  complicated  pathway-through  various  hidden  signals  that  may 
be  components  of  other  links  in  the  system.  Thus  a  link  is  associated  with  a  particular 
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set  of  dynamics-a  system-that  characterizes  how  the  input  signal  is  transformed  into  the 
output  signal.  The  fact  that  hidden  signals  may  be  shard  between  links,  however,  is  an 
important  distinction  between  signal  and  subsystem  interconnection  structures.  Note  that 
a  state  of  one  subsystem,  interconnected  with  others  in  a  subsystem  architecture  (such 
as  a  standard  feedback  interconnection  between  two  blocks),  is  never  shared  with  other 
subsystems;  the  subsystem  architecture  effectively  partitions  the  states  of  the  networked 
system.  In  contrast,  states  on  the  links  of  the  signal  structure  may,  in  fact,  be  shared  with 
those  of  other  links.  This  degree  of  abstraction  is  important  for  security  problems  because 
an  additive  perturbation  on  a  link  of  the  signal  structure  does  not  represent  the  corruption 
of  a  particular  channel,  as  it  would  in  the  subsystem  structure,  but  rather  the  idea  that  an 
attacker  infiltrated  a  particular  dependency  between  specific  manifest  variables. 

The  next  section  provides  some  background  on  Dynamical  Structure  Functions  (DSF), 
which  are  used  to  represent  the  signal  structure  of  a  system.  The  DSF  is  a  system  represen¬ 
tation  that  describes  more  structure,  in  the  logical  interconnection  sense,  than  the  transfer 
function  provides,  but  less  than  the  state  space  realization  would  reveal.  Specifically,  the 
DSF  describes  exactly  the  causal  dependencies  between  manifest  variables  without  offering 
any  indication  of  the  structure  relating  hidden  variables.  As  a  result,  although  every  state 
space  realization  specifies  a  unique  DSF,  and  every  DSF  specifies  a  unique  transfer  function, 
there  are  many  DSF  architectures  consistent  with  any  specific  transfer  function,  and  many 
state  space  realizations  consistent  with  any  specific  DSF. 

4.6.3  Background:  Dynamical  Structure  Function 

Before  developing  the  main  theorem,  we  will  present  a  concise  derivation  of  the  dynamical 
structure  function,  and  explain  its  relevance  to  the  security  of  a  networked  system.  For  a 
complete  derivation  and  results  on  different  representations  of  structure  see  [87,  88]. 

Let  us  consider  a  state-space  LTI  system 


Zi 

An 

A12 

Zi 

+ 

~B{ 

A 

A21 

A22 

Z2_ 

b2 

y  =  [C\  C2 


where  \C\  C2  has  full  row  rank.  This  system  can  be  transformed  to: 


An  A12 

A21  a22 


+ 


B2 


u 


y=  [1  0] 


(72) 


(73) 


Here  y  are  the  states  that  are  measured,  and  x  are  the  hidden  states.  Now,  taking  Laplace 
Transforms  of  the  signals  in  (73),  we  get 


sY 

An  A12 

'Y~ 

~B{ 

sX 

A21  A22 

X 

+ 

.A 

(74) 
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Solving  for  X  in  the  second  equation  of  74  gives 


A  —  {si  —  A22)  1A.2i  Y  +  ( si  —  A22)  1 B2U 
Substituting  into  the  first  equation  of  (74)  we  get, 

sY  =  WY  +  VU, 

where  W  =  An  +  Ai2(sJ  —  A22)_1A2i  and  V  =  Au^sl  —  A22)~1B2  +  Bi.  Let  D  be  a  diagonal 
matrix  with  the  diagonal  entries  of  W.  Then, 

{si  -  D)Y  —  {W  —  D)Y  +  VU. 

Now  we  can  rewrite  this  equation  as, 


Y  =  QY  +  PU,  (75) 

where 

Q  =  {si  -  D)~\W  -  D) 

and 

P={sl-  Dy'v. 

The  matrix  Q  is  a  matrix  of  transfer  functions  from  Y)  to  Y),  i  ^  j,  or  relating  each  measured 
signal  to  the  other  measured  signals.  Note  that  Q  is  zero  on  the  diagonal  and  either  zero  or 
a  strictly  proper  transfer  function  on  the  off  diagonal.  The  matrix  P  is  a  matrix  of  zeros  or 
strictly  proper  transfer  functions  from  each  input  to  each  output  without  depending  on  any 
additional  measured  states.  Together,  the  pair  (Q(s),  P(s ))  is  called  the  dynamical  structure 
function  for  system  (72). 

The  transfer  function  matrix  for  this  system  is  given  by 


G  =  {I  —  Q)-lP  =  C{sl  -  A)~lB. 


Hence,  Gtl  is  the  closed  loop  transfer  function  from  input  j  to  state  i.  In  this  section,  we 
will  also  refer  to  the  closed  loop  transfer  function  between  states.  A  transfer  function  from 
state  j  to  state  i  is  represented  by  Hl3 ,  where 

H  =  (I-Q)-\ 

Note  that  the  transfer  function  from  a  state  to  an  input  is  always  zero. 

Definition  11.  Given  a  system  12  with  signal  structure  characterized  by  the  dynamical 
structure  function  (P,Q),  a  link  {i,j)  of  the  system  corresponds  to  any  nonzero  entry  in  P 
or  Q. 

Note  that  P  gives  the  links  from  the  inputs  to  the  measured  states,  and  Q  gives  the  links 
that  represent  the  dependencies  between  the  measured  states.  The  next  section  will  introduce 
the  notion  of  vulnerability  and  characterize  vulnerable  links  in  the  system’s  architecture 
characterized  by  {P,Q). 
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5  +  2 

Figure  33:  The  system  with  the  perturbation  A.  Black  arrows  indicate  secure  links,  while 
blue  arrows  indicate  vulnerable  links. 


4.6.4  Vulnerable  Links 

In  this  work,  vulnerability  refers  to  the  destabilization  of  a  system  resulting  from  the  cor¬ 
ruption  of  a  single  link  in  its  signal  architecture.  We  begin  with  a  definition  of  a  vulnerable 
link. 


Definition  12.  Given  a  system  12  with  signal  structure  characterized  by  the  dynamical  struc¬ 
ture  function  (P,Q),  a  link  in  (P,Q)  is  called  vulnerable  if  there  exists  a  stable  perturbation 
on  the  link  that  makes  the  system  unstable. 


Example  7.  Let  us  consider  a  system  with 


-  1 

0 

0 

1  - 

p  = 

s+2 

0 

1 

s+2. 

,  and  Q  = 

1 

.s+2 

s+2 

0 

This  system  is  stable  because  the  transfer  function. 


G 


1 

s2  +  4s  +  3 


s  +  2 
1 


1 

s  +  2 


Now  let  us  add  a  perturbation  A  =  AL_  1°  the  link  Q 12  as  shown  in  Figure  33.  The  resulting 
transfer  function  is 

C  =  1  s  + 2  -1 

s(s  +  4)  [  4  s  +  2_|  ’ 

which  is  unstable.  Hence  the  link  Q 12  is  a  vulnerable  link.  Similarly,  it  can  be  shown  that 
Q 21  is  vulnerable,  although  neither  Pu  nor  P2 2  are  vulnerable. 


Condition  for  Vulnerability  Given  that  an  attacker  has  the  knowledge  of  the  dynam¬ 
ical  structure  function  representation  of  a  system,  we  will  derive  a  necessary  and  sufficient 
condition  for  a  link  to  be  vulnerable. 

Theorem  5.  Let  us  consider  a  stable  system  ( P,Q ).  There  exists  a  stable  additive  pertur¬ 
bation  A  on  a  link  from  node  i  to  node  j,  either  in  P  or  Q,  that  makes  the  system  unstable 
if  and  only  if  the  closed  loop  transfer  function  from  node  j  to  i  is  nonzero. 
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Figure  34:  System  with  the  perturbation  Ae*eJ 


Figure  35:  Necessary  and  sufficient  condition  for  stability  of  the  system  in  Figure  34 


Proof.  The  system  with  the  perturbation  A  can  be  represented  as  the  linear  fractional  trans¬ 
formation  in  Figure  34,  where  T  is  the  associated  closed  loop  transfer  function,  and  Wi ,  Wj 
represent  the  signals  at  node  i  and  j  respectively.  This  system  is  stable  if  and  only  if  the 
system  in  Figure  35  is  stable  (see  [91]).  If  T,j  =  0,  any  stable  A  does  not  affect  the  stability 
of  the  system  in  Figure  35.  Thus  the  closed  loop  system  in  Figure  34  is  stable  for  all  A. 

If  Tij  ^  0,  then  the  system  in  Figure  35  is  unstable  if  any  of  the  transfer  functions  of 


-> 


wj 

Wi 


is  unstable.  We  have, 


Wi  = 


T  A 

± ij 


PtA  A] 


Let  =  jj  and  A  =  jf,  then 


<5d 


D5 


Wj  = 


D 


DS 


D 


\NSN  <Sjv 

dj 

[  DSd  Sd 

d{ 

(76) 


For  a  polynomial  to  be  stable  it  is  necessary  that  all  its  coefficients  are  of  the  same  sign. 
In  the  case  of  the  polynomial 

R{s)  =  DSd  -  N5n,  (77) 

it  is  easy  to  see  that  a  properly  designed  A  can  zero  out  at  least  one  of  the  terms.  Thus, 
there  exists  a  A  that  destabilizes  these  transfer  functions.  □ 

Note  that  when  we  are  considering  the  vulnerability  of  the  links  in  Q,  T  =  H  =  (/— Q)-1, 
gives  the  closed  loop  transfer  functions.  Now,  we  will  present  some  implications  of  this  result. 


Corollary  2.  None  of  the  links  in  P  are  vulnerable. 
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Proof.  This  is  true  because  the  transfer  function  from  the  states  to  the  input  is  always 
zero.  □ 

Corollary  3.  If  Tl}  is  nonzero,  there  exists  a  perturbation  A  e  R  that  destabilizes  the  system 
in  Figure  35. 

Proof.  Let  A  e  M,  lt]  =  jf-  Thus,  j1  =  ADf^Nl ;  and  the  polynomial  in  (77)  becomes 
DiD  —  N(DiA  +  Ni).  We  can  see  that  at  least  one  of  the  terms  in  this  polynomial  can  be 
zeroed  out  by  choosing  appropriate  A,  making  the  polynomial  unstable.  □ 

Corollary  4.  Let  us  consider  a  stable  system, 

x  =  Ax  +  Iu ,  (78) 

V  =  lx, 

where  A  e  Mnxn  and  let  G  =  ( si  —  A)-1.  There  exists  a  perturbation  K  =  Ae,eJ,  A  e  M, 
such  that  (A  +  K)  is  not  Hurwitz,  if  and  only  if  the  transfer  function  from  input  Ui  to  output 
yj,  Gji,  is  nonzero. 

Proof.  If  the  perturbation  is  on  the  diagonal  entry  of  A,  then  it  is  easy  to  see  that  a  desta¬ 
bilizing  perturbation  always  exists  and  Gu  is  never  zero.  Let  D  =  diag(An,  A22,  •••,  Ann). 
The  dynamical  structure  function  of  the  system  is  given  by  P  =  (si  —  D)^1  and  Q  = 
(si  —  D)~l(A  —  D ).  Any  perturbation  K  =  A eiej,i  ^  j  effects  only  the  link  Qp.  Hence, 
the  perturbation  can  make  the  system  unstable  if  and  only  if  the  transfer  function  H3i  is 
nonzero.  Also,  the  diagonal  entries  of  P  are  nonzero,  and  G  =  HP.  Thus,  the  transfer 
function  H%]  is  nonzero  if  and  only  if  Gji  is  nonzero. 

□ 

Example  8.  Let  us  consider  a  system  of  the  form  78  with 

"-1  0-43" 

2  -2  0  0 

3  0  -2  -4  ' 

_  0  3-2  —  5_ 

Here  the  eigenvalue  of  A  are  a  =  {  —  1.5000  +  3. 4278j,  —1.5000  —  3. 4278j,  —6.7016,  —0.2984}. 
Hence,  the  system  is  stable.  In  this  system,  the  link  from  X4  to  x\  is  not  vulnerable  because 
G41  =  0.  Notice  that  this  example  is  not  a  trivial  example,  like  a  diagonal  or  a  triangular 
system,  since  there  are  cycles  that  contain  both  nodes  X\  and  X4. 

Corollary  5.  Let  A  G  lZnxn.  A  perturbation  on  the  (i,j)th  entry  of  A  changes  its  eigenvalues 
if  and  only  if  the  Gji  0,  where  G  =  (si  —  A)^1  is  the  transfer  function  matrix  i.e.  the 
(■ i,j )  minor  of  (si  —  A)  is  nonzero. 

Proof.  Take  the  system  from  Corollary  4.  We  can  see  that  a  perturbation  on  the  (i,j)th 
entry  has  no  effect  on  the  system  if  GtJ  =  0.  Also,  if  Gji  0,  the  perturbation  forms  a 
closed  loop  system,  such  as  the  one  given  in  Figure  35,  in  which  case  A  definitely  changes 
the  poles  of  the  system.  □ 

If  we  take  the  A  matrix  from  Example  8,  note  that  its  eigenvalues  stay  unchanged  for 
any  perturbation  on  the  (l,4)th  entry. 
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i  +  2 


Figure  36:  A  system  with  a  secure  link  in  a  cycle.  Black  arrows  represent  the  secure  links. 

Structure  and  Vulnerability  To  perform  the  vulnerability  analysis  of  a  system,  we 
assume  that  the  attacker  can  only  modify  existing  links  and  cannot  create  new  links.  With 
this  assumption,  we  can  see  that  systems  where  the  output  nodes  do  not  form  a  cycle  are 
always  secure,  because  in  such  a  case  the  nodes  can  be  permuted  to  obtain  a  triangular  Q 
matrix.  A  triangular  Q  matrix  gives  a  triangular  P,  and  by  applying  Theorem  5  we  can  see 
that  all  the  existing  links  are  secure.  Note  that  secure  links  doesn’t  always  mean  they  are 
from  a  triangular  system.  For  example,  the  link  Q14  is  secure  in  the  system  given  in  Figure 
36,  which  is  the  signal  structure  architecture  of  the  state-space  system  in  Example  8. 

Noting  that  certain  graphical  structures  result  in  secure  links  begs  the  question  of  whether 
there  are  particular  dynamics  that  contribute  to  secure  or  vulnerable  links  in  the  system’s 
architecture.  The  following  theorem  answers  this  question. 

Theorem  6.  Every  transfer  function  G  has  a  completely  secure  architecture  (. P,Q ). 

Proof.  For  any  transfer  function  G ,  note  that  (P  —  G,  Q  —  0)  is  an  admissible  Dynamical 
Structure  Function  since  G  =  (/  —  0)-1G.  From  Corollary  1,  we  see  that  none  of  the  links 
in  P  are  vulnerable,  and  since  Q  has  no  links,  the  system  is  secure.  □ 

This  result  shows  that  the  vulnerability  of  a  system  is  structure  dependent  and  not  a 
function  of  the  system  dynamics.  This  fact  highlights  one  difference  between  the  vulnerabil¬ 
ity,  which  depends  on  the  system  structure  and  not  the  dynamics,  and  the  robustness,  which 
depends  on  the  dynamics  and  not  the  system  structure. 

Measure  of  Vulnerability  Feedback  is  very  common  in  both  natural  and  engineered 
systems.  Nevertheless,  such  structures  usually  generate  vulnerable  links.  Thus,  a  measure 
of  vulnerability  is  essential  to  understand  the  security  of  the  system. 

Given  a  signal  architecture  (P,  Q)  with  associated  closed  loop  transfer  function  T,  the 
vulnerability  of  link  (i ,  j )  is  given  by 


Vji  = 


T. 

-t  7 


IJ  I  I  o° 


(79) 


which  is  the  smallest  perturbation  required  on  link  (i,j)  to  destabilize  the  system.  Since  all 
the  links  in  P  are  secure,  we  only  consider  the  links  in  Q  while  computing  the  vulnerability, 
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hence,  T  =  H.  The  vulnerability  of  the  system  is  given  by 


V  =  min  v.ji 

1 

=  nun  — — — 

(++Q  ||-7)j||oo 


(80) 

(81) 


This  measure  allows  us  to  associate  a  size  of  the  smallest  destabilizing  perturbation  with 
every  link  in  the  system  architecture.  Secure  links  thus  have  a  vulnerability  of  oo.  Note  that 
V,  the  system  vulnerability,  is  different  than  (and  not  smaller  than)  the  size  of  the  smallest 
destabilizing  perturbation  for  the  system,  since  link  perturbations  are  restricted  to  act  on  a 
single  link  only. 


4.6.5  Numerical  Example 


Let  us  consider  a  system  with  the  architecture  given  in  Figure  37a  where, 


P  = 


r  i 


S+1 

0 

0 


0 

1 

s+1 

0 


0  ' 

0 

1 

s+1. 


and 


Q 


o 

i 

s+2 

0 


0 

0 

1 

s+3 


s+ll 


0 

0 


The  transfer  function  matrix  for  the  system  is  given  by 


G 


~s3+6s2+11s+6  s+2  s2+5s+6 

d(s)  d(s)  d(s) 

s2+4s+3  s3+6s2+11s+6  s+3 

d(s)  d(s)  d(s) 

s+1  s2+3s+2  s3+6s2+11s+6 

d(s)  d(s)  d(s) 


where  d(s)  =  s4  +  7s3  +  17s2  +  16s  +  5.  The  gain  of  this  system  is  given  by  l+Hoo  =  2.4. 
Hence,  by  small  gain  theorem,  a  perturbation  of  gain  larger  than  0.42  could  destabilize  the 
system. 

Let  H  =  (7  — Q)_1  represent  the  transfer  function  between  the  measured  states  yt.  Since 
the  links  in  P  are  not  vulnerable,  we  consider  the  perturbations  on  the  links  in  Q  which  are 
the  links  (2/1?  2/2) 5  (' 2/2, 2/3 ),  and  (2/3, 7/i) -  To  compute  the  vulnerability  of  these  links  we  need 
the  following  transfer  functions: 


s  +  2 

s3  +  6s2  +  11s  +  5 
s  +  3 

s3  +  6s2  +  11s  +  5 
s  +  1 

s3  +  6s2  +  11s  +  5 

For  this  system  v12  =  Ai2  =  2.5,  V23  =  A23  =  |,  and  V31  =  A31  =  5.  Hence,  the  smallest 
perturbation  on  a  single  link  that  can  destabilize  this  system  must  have  a  gain  of  1.67. 


H 12  = 

7723  = 
773i  = 
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Figure  37:  Vulnerable  and  secure  architectures  for  the  same  transfer  function.  Black  links 
are  secure,  vulnerable  links  are  colored  blue,  yellow,  and  red  in  the  increasing  order  of  their 
vulnerability. 


This  system  can  also  be  implemented  as  shown  in  Figure  37b,  where  P  —  G.  This  is 
one  of  the  secure  implementations  of  the  system  in  Figure  37a.  From  this  example  we  thus 
observe  the  following: 

•  The  same  transfer  function  can  exhibit  both  vulnerable  and  secure  architectures, 

•  System  robustness,  characterized  by  the  size  of  the  smallest  destabilizing  perturbation 
(0.42  in  this  example),  is  not  equivalent  to  system  vulnerability,  characterized  by  the 
size  of  the  smallest  destabilizing  perturbation  on  a  single  link  (about  1.67  in  this 
example) , 

•  Only  links  in  Q  can  be  vulnerable. 

4.6.6  Summary 

This  section  explored  the  notion  of  a  vulnerable  link  in  a  network  of  controlled  linear  dynam¬ 
ical  systems.  The  architecture  of  the  system  was  characterized  by  its  Dynamical  Structure 
Function,  representing  the  logical  interconnection  structure  of  the  system.  Vunerability  was 
then  defined  as  the  size  of  the  smallest  destabilizing  perturbation  acting  on  a  single  link. 
The  main  results  of  the  section  provided  necessary  and  sufficient  conditions  for  the  vulner¬ 
ability  of  a  link  and  then  demonstrated  that  any  transfer  function  has  a  completely  secure 
architecture.  This  result  highlights  the  idea  that  while  robustness  is  a  property  of  a  system’s 
dynamics,  security  (in  the  sense  discussed  here)  is  a  property  of  its  signal  architecture.  Fu¬ 
ture  work  will  explore  parallel  notions  for  subsystem  interconnection  structure  of  networked 
systems. 
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5  Conclusions 


This  work  considered  a  new  model  for  complex  network  environments  and  explored  its  ap¬ 
plicability  to  wireless  mesh  networks  and  bio-chemical  reaction  networks.  The  new  model 
characterizes  structure  of  a  networked  system  in  a  novel  way  that  does  not  attempt  to  parti¬ 
tion  every  state  variable  of  the  system  into  well  defined  subsystems,  resulting  in  a  significant 
reduction  in  the  information  load  necessary  to  reconstruct  the  (new)  network  structure  from 
experimental  data.  Instead,  the  new  model  characterizes  the  structure  of  the  network  as 
the  dependency  among  manifest  variables  (i.e.  those  variables  that  are  visible  to  the  outside 
world),  and  this  structure  can  be  discovered  with  O(n)  experiments,  where  n  is  the  number 
of  manifest  variables  (observed  nodes).  Moreover,  algorithms  for  network  reconstruction 
were  developed  that  apply  to  nonlinear  systems  near  equilibrium  with  noisy  measurements. 
These  “robust”  reconstruction  methods  have  been  tested  in-silico,  and  experiments  to  vali¬ 
date  them  on  both  wireless  mesh  and  bio-chemical  reaction  networks  remain  a  goal  for  future 
research. 

Although  the  new  network  model  can  be  applied  directly  to  existing  models  of  chemical 
reaction  kinetics,  application  to  wireless  mesh  networks  is  less  obvious.  As  a  result,  significant 
efforts  were  made  to  model  and  understand  the  behavior  of  wireless  mesh  networks  leveraging 
the  team’s  varied  expertise  through  the  “network  design  cycle.”  Following  this  process,  we 
introduced  a  new  partial  interference  model  and  a  new  first  principles  model  for  wireless  mesh 
networks,  used  them  to  design  novel  rate  control  protocols,  and  then  verified  the  improved 
performance  of  these  protocols  through  both  simulation  and  experimental  testing.  Although 
these  models  are  not  yet  exactly  the  type  of  model  needed  for  network  reconstruction,  they  lay 
the  appropriate  foundation  for  developing  such  a  model.  Our  next  steps  will  be  to  extend  this 
work  to  a  model  appropriate  for  network  reconstruction  in  wireless  mesh  networks,  and  then 
use  this  model  to  experimentally  validate  our  reconstruction  methods  for  these  networked 
systems. 

The  new  network  reconstruction  model  also  enables  a  rigorous  analysis  of  network  vul¬ 
nerabilities.  Preliminary  work  has  demonstrated  that  certain  network  links  can  exhibit 
system-wide  vulnerability,  and  corresponding  analysis  has  demonstrated  that  every  system 
has  a  completely  secure  architecture.  In  practice  this  architecture  may  be  unrealizable,  since 
feedback  is  often  an  essential  part  of  networked  systems,  but  these  early  results  point  to 
the  power  of  the  new  models  in  enabling  a  rigorous  vulnerability  analysis  by  exploring  the 
robustness  of  the  closed-loop  system.  Future  work  will  build  on  these  results  and  continue  to 
explore  security  implications  of  network  attacks.  Ultimately,  for  example,  we  would  like  to 
be  able  to  exploit  our  network  reconstruction  technologies  to  detect  intruders  in  our  wireless 
mesh  networks  and  understand  how  to  deploy  specific  topologies  that  minimizes  security 
risks  of  link-wise  attacks. 
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